Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MappingGeneral.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: MappingGeneral.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: Mapping of the type X = X(x,y), Y = Y(x,y)
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 #include <MultiRegions/ExpList.h>
38 #include <iomanip>
39 
40 namespace Nektar
41 {
42 namespace GlobalMapping
43 {
44 
45 std::string MappingGeneral::className =
47  MappingGeneral::create, "X = X(x,y,z), Y = Y(x,y,z), Z=Z(x,y,z)");
48 
49 /**
50  * @class MappingGeneral
51  * This class implements the most general mapping, defined by the transformation
52  * \f[ \bar{x} = \bar{x}(x,y,z) \f]
53  * \f[ \bar{y} = \bar{y}(x,y,z) \f]
54  * \f[ \bar{z} = \bar{z}(x,y,z) \f]
55  * where \f$(\bar{x},\bar{y},\bar{z})\f$ are the Cartesian (physical)
56  * coordinates and \f$(x,y,z)\f$ are the transformed (computational)
57  * coordinates.
58  */
62  : Mapping(pSession, pFields)
63 {
64 }
65 
66 
67 /**
68  *
69  */
72  const TiXmlElement *pMapping)
73 {
74  Mapping::v_InitObject(pFields, pMapping);
75 
76  m_constantJacobian = false;
77 
79  "General Mapping needs at least 2 velocity components.");
80 }
81 
83  const Array<OneD, Array<OneD, NekDouble> > &inarray,
84  Array<OneD, Array<OneD, NekDouble> > &outarray)
85 {
86  int physTot = m_fields[0]->GetTotPoints();
87  int nvel = m_nConvectiveFields;
88 
89  for (int i=0; i<nvel; i++)
90  {
91  outarray[i] = Array<OneD, NekDouble> (physTot, 0.0);
92  for (int j=0; j<nvel; j++)
93  {
94  Vmath::Vvtvp(physTot, inarray[j], 1,
95  m_deriv[i*nvel+j], 1,
96  outarray[i], 1,
97  outarray[i], 1);
98  }
99 
100  }
101 }
102 
104  const Array<OneD, Array<OneD, NekDouble> > &inarray,
105  Array<OneD, Array<OneD, NekDouble> > &outarray)
106 {
107  int physTot = m_fields[0]->GetTotPoints();
108  int nvel = m_nConvectiveFields;
109 
110  for (int i=0; i<nvel; i++)
111  {
112  outarray[i] = Array<OneD, NekDouble> (physTot, 0.0);
113  for (int j=0; j<nvel; j++)
114  {
115  Vmath::Vvtvp(physTot, inarray[j], 1,
116  m_invDeriv[i*nvel+j], 1,
117  outarray[i], 1,
118  outarray[i], 1);
119  }
120 
121  }
122 }
123 
125  const Array<OneD, Array<OneD, NekDouble> > &inarray,
126  Array<OneD, Array<OneD, NekDouble> > &outarray)
127 {
128  int physTot = m_fields[0]->GetTotPoints();
129  int nvel = m_nConvectiveFields;
130 
131  for (int i=0; i<nvel; i++)
132  {
133  outarray[i] = Array<OneD, NekDouble> (physTot, 0.0);
134  for (int j=0; j<nvel; j++)
135  {
136  Vmath::Vvtvp(physTot, inarray[j], 1,
137  m_invDeriv[j*nvel+i], 1,
138  outarray[i], 1,
139  outarray[i], 1);
140  }
141 
142  }
143 }
144 
146  const Array<OneD, Array<OneD, NekDouble> > &inarray,
147  Array<OneD, Array<OneD, NekDouble> > &outarray)
148 {
149  int physTot = m_fields[0]->GetTotPoints();
150  int nvel = m_nConvectiveFields;
151 
152  for (int i=0; i<nvel; i++)
153  {
154  outarray[i] = Array<OneD, NekDouble> (physTot, 0.0);
155  for (int j=0; j<nvel; j++)
156  {
157  Vmath::Vvtvp(physTot, inarray[j], 1,
158  m_deriv[j*nvel+i], 1,
159  outarray[i], 1,
160  outarray[i], 1);
161  }
162 
163  }
164 }
165 
167  Array<OneD, NekDouble> &outarray)
168 {
169  int physTot = m_fields[0]->GetTotPoints();
170  Vmath::Vcopy(physTot, m_jac, 1, outarray, 1);
171 }
172 
174  Array<OneD, Array<OneD, NekDouble> > &outarray)
175 {
176  int physTot = m_fields[0]->GetTotPoints();
177  int nvel = m_nConvectiveFields;
178 
179  for (int i=0; i<nvel; i++)
180  {
181  for (int j=0; j<nvel; j++)
182  {
183  outarray[i*nvel+j] = Array<OneD, NekDouble> (physTot, 0.0);
184  Vmath::Vcopy(physTot, m_metricTensor[i*nvel+j], 1,
185  outarray[i*nvel+j], 1);
186  }
187 
188  }
189 }
190 
192  Array<OneD, Array<OneD, NekDouble> > &outarray)
193 {
194  int physTot = m_fields[0]->GetTotPoints();
195  int nvel = m_nConvectiveFields;
196 
197  for (int i=0; i<nvel; i++)
198  {
199  for (int j=0; j<nvel; j++)
200  {
201  outarray[i*nvel+j] = Array<OneD, NekDouble> (physTot, 0.0);
202  Vmath::Vcopy(physTot, m_invMetricTensor[i*nvel+j], 1,
203  outarray[i*nvel+j], 1);
204  }
205 
206  }
207 }
208 
210  const Array<OneD, Array<OneD, NekDouble> > &inarray,
211  Array<OneD, Array<OneD, NekDouble> > &outarray)
212 {
213  int physTot = m_fields[0]->GetTotPoints();
214  int nvel = m_nConvectiveFields;
215 
216  // Calculate {i,jk} u^j
217  for (int i = 0; i <nvel; i++)
218  {
219  for (int k = 0; k < nvel; k++)
220  {
221  outarray[i*nvel+k] = Array<OneD, NekDouble>(physTot,0.0);
222  for (int j = 0; j < nvel; j++)
223  {
224  Vmath::Vvtvp(physTot, inarray[j], 1,
225  m_Christoffel[i*nvel*nvel+j*nvel+k], 1,
226  outarray[i*nvel+k], 1,
227  outarray[i*nvel+k], 1);
228  }
229  }
230  }
231 
232 }
233 
235  const Array<OneD, Array<OneD, NekDouble> > &inarray,
236  Array<OneD, Array<OneD, NekDouble> > &outarray)
237 {
238  int physTot = m_fields[0]->GetTotPoints();
239  int nvel = m_nConvectiveFields;
240 
241  // Calculate {i,jk} u_i
242  for (int j = 0; j < nvel; j++)
243  {
244  for (int k = 0; k < nvel; k++)
245  {
246  outarray[j*nvel+k] = Array<OneD, NekDouble>(physTot,0.0);
247  for (int i = 0; i <nvel; i++)
248  {
249  Vmath::Vvtvp(physTot, inarray[i], 1,
250  m_Christoffel[i*nvel*nvel+j*nvel+k], 1,
251  outarray[j*nvel+k], 1,
252  outarray[j*nvel+k], 1);
253  }
254  }
255  }
256 }
257 
259 {
262 }
263 
265 {
266  int physTot = m_fields[0]->GetTotPoints();
267  int nvel = m_nConvectiveFields;
268 
269  // Set wavespace to false and store current value
270  bool wavespace = m_fields[0]->GetWaveSpace();
271  m_fields[0]->SetWaveSpace(false);
272 
273  // Allocate memory
278  for (int i = 0; i < m_metricTensor.num_elements(); i++)
279  {
280  m_metricTensor[i] = Array<OneD, NekDouble>(physTot, 0.0);
281  m_invMetricTensor[i] = Array<OneD, NekDouble>(physTot, 0.0);
282  m_deriv[i] = Array<OneD, NekDouble>(physTot, 0.0);
283  m_invDeriv[i] = Array<OneD, NekDouble>(physTot, 0.0);
284  }
285  m_jac = Array<OneD, NekDouble>(physTot, 0.0);
286 
287  // First, calculate derivatives of the mapping -> dX^i/dx^j = c^i_j
288  for( int i = 0; i<nvel; i++)
289  {
290  for( int j = 0; j<nvel; j++)
291  {
292  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[j],
293  m_coords[i],
294  m_deriv[i*nvel+j]);
295  }
296  }
297  // In Homogeneous case, m_deriv(2,2) needs to be set to 1
298  // because differentiation in wavespace is not valid for non-periodic field
299  if (m_fields[0]->GetExpType() == MultiRegions::e3DH1D)
300  {
301  Vmath::Fill(physTot, 1.0, m_deriv[2*nvel+2], 1);
302  }
303 
304  // Now calculate the metric tensor --> g_ij = sum_k { c^k_i c^k_j }
305  for( int i = 0; i<nvel; i++)
306  {
307  for( int j = 0; j<nvel; j++)
308  {
309  for( int k = 0; k<nvel; k++)
310  {
311  Vmath::Vvtvp(physTot, m_deriv[k*nvel+i], 1,
312  m_deriv[k*nvel+j], 1,
313  m_metricTensor[i*nvel+j], 1,
314  m_metricTensor[i*nvel+j], 1);
315  }
316  }
317  }
318 
319  // Put the adjoint of g in m_invMetricTensor
320  switch (nvel)
321  {
322  case 1:
323  Vmath::Fill (physTot, 1.0, m_invMetricTensor[0], 1);
324  break;
325  case 2:
326  Vmath::Vcopy(physTot, m_metricTensor[1*nvel+1], 1,
327  m_invMetricTensor[0*nvel+0], 1);
328  Vmath::Smul (physTot, -1.0, m_metricTensor[0*nvel+1], 1,
329  m_invMetricTensor[1*nvel+0], 1);
330  Vmath::Smul (physTot, -1.0, m_metricTensor[1*nvel+0], 1,
331  m_invMetricTensor[0*nvel+1], 1);
332  Vmath::Vcopy(physTot, m_metricTensor[0*nvel+0], 1,
333  m_invMetricTensor[1*nvel+1], 1);
334  break;
335  case 3:
336  {
337  int a, b, c, d, e, i, j;
338 
339  // Compute g^{ij} by computing Cofactors(g_ij)^T
340  for (i = 0; i < nvel; ++i)
341  {
342  for (j = 0; j < nvel; ++j)
343  {
344  a = ((i+1)%nvel) * nvel + ((j+1)%nvel);
345  b = ((i+1)%nvel) * nvel + ((j+2)%nvel);
346  c = ((i+2)%nvel) * nvel + ((j+1)%nvel);
347  d = ((i+2)%nvel) * nvel + ((j+2)%nvel);
348  e = i*nvel + j;
349  // a*d - b*c
350  Vmath::Vmul(physTot, m_metricTensor[b], 1,
351  m_metricTensor[c], 1,
352  m_invMetricTensor[e], 1);
353  Vmath::Vvtvm(physTot, m_metricTensor[a], 1,
354  m_metricTensor[d], 1,
355  m_invMetricTensor[e], 1,
356  m_invMetricTensor[e], 1);
357  }
358  }
359  break;
360  }
361  }
362 
363  // Compute g = det(g_{ij}) (= Jacobian squared) and store
364  // temporarily in m_jac.
365  for (int i = 0; i < nvel; ++i)
366  {
367  Vmath::Vvtvp(physTot, m_metricTensor[i], 1,
368  m_invMetricTensor[i*nvel], 1,
369  m_jac, 1, m_jac, 1);
370  }
371 
372  // Calculate g^ij (the inverse of g_ij) by dividing by jac
373  for (int i = 0; i < nvel*nvel; ++i)
374  {
375  Vmath::Vdiv(physTot, m_invMetricTensor[i], 1, m_jac, 1,
376  m_invMetricTensor[i], 1);
377  }
378 
379  // Compute the Jacobian = sqrt(g)
380  Vmath::Vsqrt(physTot, m_jac, 1, m_jac, 1);
381 
382  // Calculate the derivatives of the inverse transformation
383  // c'^j_i = dx^j/dX^i = sum_k {g^jk c^i_k}
384  for (int i = 0; i < nvel; ++i)
385  {
386  for (int j = 0; j < nvel; ++j)
387  {
388  for (int k = 0; k < nvel; ++k)
389  {
390  Vmath::Vvtvp(physTot, m_deriv[i*nvel+k], 1,
391  m_invMetricTensor[j*nvel+k], 1,
392  m_invDeriv[i*nvel+j], 1,
393  m_invDeriv[i*nvel+j], 1);
394  }
395  }
396  }
397 
398  // Restore value of wavespace
399  m_fields[0]->SetWaveSpace(wavespace);
400 }
401 
403 {
404  int physTot = m_fields[0]->GetTotPoints();
405  int nvel = m_nConvectiveFields;
406 
407  Array<OneD, Array<OneD, NekDouble> > gradG(nvel*nvel*nvel);
408  Array<OneD, Array<OneD, NekDouble> > tmp(nvel*nvel*nvel);
410  // Allocate memory
411  for (int i = 0; i < gradG.num_elements(); i++)
412  {
413  gradG[i] = Array<OneD, NekDouble>(physTot, 0.0);
414  tmp[i] = Array<OneD, NekDouble>(physTot, 0.0);
415  m_Christoffel[i] = Array<OneD, NekDouble>(physTot, 0.0);
416  }
417 
418  // Set wavespace to false and store current value
419  bool waveSpace = m_fields[0]->GetWaveSpace();
420  m_fields[0]->SetWaveSpace(false);
421 
422  //Calculate gradients of g_ij
423  for (int i = 0; i <nvel; i++)
424  {
425  for(int j=0; j<nvel; j++)
426  {
427  for(int k=0; k<nvel; k++)
428  {
429  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[k],
430  m_metricTensor[i*nvel+j],
431  gradG[i*nvel*nvel + j*nvel + k]);
432  }
433  }
434  }
435 
436  // Calculate tmp[p,j,k] = 1/2( gradG[pj,k]+ gradG[pk,j]-gradG[jk,p])
437  for (int p = 0; p <nvel; p++)
438  {
439  for (int j = 0; j < nvel; j++)
440  {
441  for (int k = 0; k < nvel; k++)
442  {
443  Vmath::Vadd(physTot, gradG[p*nvel*nvel + j*nvel + k], 1,
444  gradG[p*nvel*nvel + k*nvel + j], 1,
445  tmp[p*nvel*nvel + j*nvel + k], 1);
446  Vmath::Vsub(physTot, tmp[p*nvel*nvel + j*nvel + k], 1,
447  gradG[j*nvel*nvel + k*nvel + p], 1,
448  tmp[p*nvel*nvel + j*nvel + k], 1);
449  Vmath::Smul(physTot, 0.5, tmp[p*nvel*nvel + j*nvel + k], 1,
450  tmp[p*nvel*nvel + j*nvel + k], 1);
451  }
452  }
453  }
454 
455  // Calculate Christoffel symbols = g^ip tmp[p,j,k]
456  for (int i = 0; i <nvel; i++)
457  {
458  for (int j = 0; j < nvel; j++)
459  {
460  for (int k = 0; k < nvel; k++)
461  {
462  for (int p = 0; p < nvel; p++)
463  {
464  Vmath::Vvtvp(physTot, m_invMetricTensor[i*nvel+p], 1,
465  tmp[p*nvel*nvel + j*nvel + k], 1,
466  m_Christoffel[i*nvel*nvel+j*nvel+k], 1,
467  m_Christoffel[i*nvel*nvel+j*nvel+k], 1);
468  }
469  }
470  }
471  }
472  // Restore wavespace
473  m_fields[0]->SetWaveSpace(waveSpace);
474 
475 }
476 
477 }
478 }
Array< OneD, Array< OneD, NekDouble > > m_invMetricTensor
Array< OneD, Array< OneD, NekDouble > > m_coords
Array with the Cartesian coordinates.
Definition: Mapping.h:411
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
virtual GLOBAL_MAPPING_EXPORT void v_CovarFromCartesian(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
MappingFactory & GetMappingFactory()
Declaration of the mapping factory singleton.
Definition: Mapping.cpp:49
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:394
virtual GLOBAL_MAPPING_EXPORT void v_ContravarFromCartesian(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Array< OneD, Array< OneD, NekDouble > > m_invDeriv
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
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
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
virtual GLOBAL_MAPPING_EXPORT void v_ApplyChristoffelCovar(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
int m_nConvectiveFields
Number of velocity components.
Definition: Mapping.h:417
virtual GLOBAL_MAPPING_EXPORT void v_GetInvMetricTensor(Array< OneD, Array< OneD, NekDouble > > &outarray)
virtual GLOBAL_MAPPING_EXPORT void v_InitObject(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const TiXmlElement *pMapping)
virtual GLOBAL_MAPPING_EXPORT void v_ContravarToCartesian(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
virtual GLOBAL_MAPPING_EXPORT void v_ApplyChristoffelContravar(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Definition: Mapping.h:409
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
Array< OneD, Array< OneD, NekDouble > > m_deriv
static std::string className
Name of the class.
Array< OneD, NekDouble > m_jac
virtual GLOBAL_MAPPING_EXPORT void v_InitObject(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const TiXmlElement *pMapping)
Definition: Mapping.cpp:100
bool m_constantJacobian
Flag defining if the Jacobian is constant.
Definition: Mapping.h:427
virtual GLOBAL_MAPPING_EXPORT void v_GetMetricTensor(Array< OneD, Array< OneD, NekDouble > > &outarray)
virtual GLOBAL_MAPPING_EXPORT void v_CovarToCartesian(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Array< OneD, Array< OneD, NekDouble > > m_Christoffel
Base class for mapping to be applied to the coordinate system.
Definition: Mapping.h:69
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
static GLOBAL_MAPPING_EXPORT MappingSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const TiXmlElement *pMapping)
Creates an instance of this class.
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:86
MappingGeneral(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
void Vvtvm(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)
vvtvm (vector times vector plus vector): z = w*x - y
Definition: Vmath.cpp:451
virtual GLOBAL_MAPPING_EXPORT void v_GetJacobian(Array< OneD, NekDouble > &outarray)
Array< OneD, Array< OneD, NekDouble > > m_metricTensor
virtual GLOBAL_MAPPING_EXPORT void v_UpdateGeomInfo()
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
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
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215