Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GeomFactors.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: GeomFactors.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: Geometric factors base class.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <SpatialDomains/GeomFactors.h>
37 #include <LibUtilities/Foundations/Interp.h>
38 
39 namespace Nektar
40 {
41  namespace SpatialDomains
42  {
43  /**
44  * @class GeomFactors
45  *
46  * This class stores the various geometric factors associated with a
47  * specific element, necessary for fundamental integration and
48  * differentiation operations as well as edge and surface normals.
49  *
50  * Initially, these algorithms are provided with a mapping from the
51  * reference region element to the physical element. Practically, this
52  * is represented using a corresponding reference region element for
53  * each coordinate component. Note that for straight-sided elements,
54  * these elements will be of linear order. Curved elements are
55  * represented using higher-order coordinate mappings. This geometric
56  * order is in contrast to the order of the spectral/hp expansion order
57  * on the element.
58  *
59  * For application of the chain rule during differentiation we require
60  * the partial derivatives \f[\frac{\partial \xi_i}{\partial \chi_j}\f]
61  * evaluated at the physical points of the expansion basis. We also
62  * construct the inverse metric tensor \f$g^{ij}\f$ which, in the case
63  * of a domain embedded in a higher-dimensional space, supports the
64  * conversion of covariant quantities to contravariant quantities.
65  * When the expansion dimension is equal to the coordinate dimension the
66  * Jacobian of the mapping \f$\chi_j\f$ is a square matrix and
67  * consequently the required terms are the entries of the inverse of the
68  * Jacobian. However, in general this is not the case, so we therefore
69  * implement the construction of these terms following the derivation
70  * in Cantwell, et. al. \cite CaYaKiPeSh13. Given the coordinate maps
71  * \f$\chi_i\f$, this comprises of five steps
72  * -# Compute the terms of the Jacobian
73  * \f$\frac{\partial \chi_i}{\partial \xi_j}\f$.
74  * -# Compute the metric tensor
75  * \f$g_{ij}=\mathbf{t}_{(i)}\cdot\mathbf{t}_{(j)}\f$.
76  * -# Compute the square of the Jacobian determinant
77  * \f$g=|\mathbf{g}|\f$.
78  * -# Compute the inverse metric tensor \f$g^{ij}\f$.
79  * -# Compute the terms \f$\frac{\partial \xi_i}{\partial \chi_j}\f$.
80  */
81 
82  /**
83  * @param gtype Specified whether the geometry is regular or
84  * deformed.
85  * @param coordim Specifies the dimension of the coordinate
86  * system.
87  * @param Coords Coordinate maps of the element.
88  */
89  GeomFactors::GeomFactors(
90  const GeomType gtype,
91  const int coordim,
92  const StdRegions::StdExpansionSharedPtr &xmap,
93  const Array<OneD, Array<OneD, NekDouble> > &coords) :
94  m_type(gtype),
95  m_expDim(xmap->GetShapeDimension()),
96  m_coordDim(coordim),
97  m_valid(true),
98  m_xmap(xmap),
99  m_coords(coords)
100  {
101  CheckIfValid();
102  }
103 
104 
105  /**
106  * @param S An instance of a GeomFactors class from which
107  * to construct a new instance.
108  */
109  GeomFactors::GeomFactors(const GeomFactors &S) :
110  m_type(S.m_type),
111  m_expDim(S.m_expDim),
112  m_coordDim(S.m_coordDim),
113  m_valid(S.m_valid),
114  m_xmap(S.m_xmap),
115  m_coords(S.m_coords)
116  {
117  }
118 
119 
120  /**
121  *
122  */
123  GeomFactors::~GeomFactors()
124  {
125  }
126 
127 
128  /**
129  * Member data equivalence is tested in the following order: shape type,
130  * expansion dimension, coordinate dimension and coordinates.
131  */
132  bool operator==(const GeomFactors &lhs, const GeomFactors &rhs)
133  {
134  if(!(lhs.m_type == rhs.m_type))
135  {
136  return false;
137  }
138 
139  if(!(lhs.m_expDim == rhs.m_expDim))
140  {
141  return false;
142  }
143 
144  if(!(lhs.m_coordDim == rhs.m_coordDim))
145  {
146  return false;
147  }
148 
149  const Array<OneD, const NekDouble> jac_lhs =
150  lhs.ComputeJac(lhs.m_xmap->GetPointsKeys());
151  const Array<OneD, const NekDouble> jac_rhs =
152  rhs.ComputeJac(rhs.m_xmap->GetPointsKeys());
153  if(!(jac_lhs == jac_rhs))
154  {
155  return false;
156  }
157 
158  return true;
159  }
160 
161 
162  /**
163  * Derivatives are computed at the geometry point distributions and
164  * interpolated to the target point distributions.
165  *
166  * @param tpoints Target point distributions.
167  * @returns Derivative of coordinate map evaluated at
168  * target point distributions.
169  */
170  DerivStorage GeomFactors::ComputeDeriv(
171  const LibUtilities::PointsKeyVector &keyTgt) const
172  {
173  ASSERTL1(keyTgt.size() == m_expDim,
174  "Dimension of target point distribution does not match "
175  "expansion dimension.");
176 
177  int i = 0, j = 0;
178  int nqtot_map = 1;
179  int nqtot_tbasis = 1;
180  DerivStorage deriv = DerivStorage(m_expDim);
181  DerivStorage d_map = DerivStorage(m_expDim);
182  LibUtilities::PointsKeyVector map_points(m_expDim);
183 
184  // Allocate storage and compute number of points
185  for (i = 0; i < m_expDim; ++i)
186  {
187  map_points[i] = m_xmap->GetBasis(i)->GetPointsKey();
188  nqtot_map *= map_points[i].GetNumPoints();
189  nqtot_tbasis *= keyTgt[i].GetNumPoints();
190  deriv[i] = Array<OneD, Array<OneD,NekDouble> >(m_coordDim);
191  d_map[i] = Array<OneD, Array<OneD,NekDouble> >(m_coordDim);
192  }
193 
194  // Calculate local derivatives
195  for(i = 0; i < m_coordDim; ++i)
196  {
197  Array<OneD, NekDouble> tmp(nqtot_map);
198  // Transform from coefficient space to physical space
199  m_xmap->BwdTrans(m_coords[i], tmp);
200 
201  // Allocate storage and take the derivative (calculated at the
202  // points as specified in 'Coords')
203  for (j = 0; j < m_expDim; ++j)
204  {
205  d_map[j][i] = Array<OneD,NekDouble>(nqtot_map);
206  deriv[j][i] = Array<OneD,NekDouble>(nqtot_tbasis);
207  m_xmap->StdPhysDeriv(j, tmp, d_map[j][i]);
208  }
209  }
210 
211  for (i = 0; i < m_coordDim; ++i)
212  {
213  // Interpolate the derivatives:
214  // - from the points as defined in the mapping ('Coords')
215  // - to the points we at which we want to know the metrics
216  // ('tbasis')
217  bool same = true;
218  for (j = 0; j < m_expDim; ++j)
219  {
220  same = same && (map_points[j] == keyTgt[j]);
221  }
222  if( same )
223  {
224  for (j = 0; j < m_expDim; ++j)
225  {
226  deriv[j][i] = d_map[j][i];
227  }
228  }
229  else
230  {
231  for (j = 0; j < m_expDim; ++j)
232  {
233  Interp(map_points, d_map[j][i], keyTgt, deriv[j][i]);
234  }
235  }
236  }
237 
238  return deriv;
239  }
240 
241 
242  /**
243  * This routine returns an array of values specifying the Jacobian
244  * of the mapping at quadrature points in the element. The array
245  * is either of size 1 in the case of elements having #GeomType
246  * #eRegular, or of size equal to the number of quadrature points for
247  * #eDeformed elements.
248  *
249  * @returns Array containing the Jacobian of the coordinate
250  * mapping at the quadrature points of the element.
251  * @see GeomType
252  */
253  Array<OneD, NekDouble> GeomFactors::ComputeJac(
254  const LibUtilities::PointsKeyVector &keyTgt) const
255  {
256  ASSERTL1(keyTgt.size() == m_expDim,
257  "Dimension of target point distribution does not match "
258  "expansion dimension.");
259 
260  int i = 0, j = 0, k = 0, l = 0;
261  int ptsTgt = 1;
262 
263  if (m_type == eDeformed)
264  {
265  // Allocate storage and compute number of points
266  for (i = 0; i < m_expDim; ++i)
267  {
268  ptsTgt *= keyTgt[i].GetNumPoints();
269  }
270  }
271 
272  // Get derivative at geometry points
273  DerivStorage deriv = ComputeDeriv(keyTgt);
274 
275  Array<TwoD, NekDouble> tmp (m_expDim*m_expDim, ptsTgt, 0.0);
276  Array<TwoD, NekDouble> gmat(m_expDim*m_expDim, ptsTgt, 0.0);
277  Array<OneD, NekDouble> jac (ptsTgt, 0.0);
278 
279  // Compute g_{ij} as t_i \cdot t_j and store in tmp
280  for (i = 0, l = 0; i < m_expDim; ++i)
281  {
282  for (j = 0; j < m_expDim; ++j, ++l)
283  {
284  for (k = 0; k < m_coordDim; ++k)
285  {
286  Vmath::Vvtvp(ptsTgt, &deriv[i][k][0], 1,
287  &deriv[j][k][0], 1,
288  &tmp[l][0], 1,
289  &tmp[l][0], 1);
290  }
291  }
292  }
293 
294  Adjoint(tmp, gmat);
295 
296  // Compute g = det(g_{ij}) (= Jacobian squared) and store
297  // temporarily in m_jac.
298  for (i = 0; i < m_expDim; ++i)
299  {
300  Vmath::Vvtvp(ptsTgt, &tmp[i][0], 1, &gmat[i*m_expDim][0], 1,
301  &jac[0], 1, &jac[0], 1);
302  }
303 
304  // Compute the Jacobian = sqrt(g)
305  Vmath::Vsqrt(ptsTgt, &jac[0], 1, &jac[0], 1);
306 
307  return jac;
308  }
309 
310 
311  /**
312  * This routine returns a two-dimensional array of values specifying
313  * the inverse metric terms associated with the coordinate mapping of
314  * the corresponding reference region to the physical element. These
315  * terms correspond to the \f$g^{ij}\f$ terms in \cite CaYaKiPeSh13 and,
316  * in the case of an embedded manifold, map covariant quantities to
317  * contravariant quantities. The leading index of the array is the index
318  * of the term in the tensor numbered as
319  * \f[\left(\begin{array}{ccc}
320  * 0 & 1 & 2 \\
321  * 1 & 3 & 4 \\
322  * 2 & 4 & 5
323  * \end{array}\right)\f].
324  * The second dimension is either of size 1 in the case of elements
325  * having #GeomType #eRegular, or of size equal to the number of
326  * quadrature points for #eDeformed elements.
327  *
328  * @see [Wikipedia "Covariance and Contravariance of Vectors"]
329  * (http://en.wikipedia.org/wiki/Covariance_and_contravariance_of_vectors)
330  * @returns Two-dimensional array containing the inverse
331  * metric tensor of the coordinate mapping.
332  */
333  Array<TwoD, NekDouble> GeomFactors::ComputeGmat(
334  const LibUtilities::PointsKeyVector &keyTgt) const
335  {
336  ASSERTL1(keyTgt.size() == m_expDim,
337  "Dimension of target point distribution does not match "
338  "expansion dimension.");
339 
340  int i = 0, j = 0, k = 0, l = 0;
341  int ptsTgt = 1;
342 
343  if (m_type == eDeformed)
344  {
345  // Allocate storage and compute number of points
346  for (i = 0; i < m_expDim; ++i)
347  {
348  ptsTgt *= keyTgt[i].GetNumPoints();
349  }
350  }
351 
352  // Get derivative at geometry points
353  DerivStorage deriv = ComputeDeriv(keyTgt);
354 
355  Array<TwoD, NekDouble> tmp (m_expDim*m_expDim, ptsTgt, 0.0);
356  Array<TwoD, NekDouble> gmat(m_expDim*m_expDim, ptsTgt, 0.0);
357  Array<OneD, NekDouble> jac (ptsTgt, 0.0);
358 
359  // Compute g_{ij} as t_i \cdot t_j and store in tmp
360  for (i = 0, l = 0; i < m_expDim; ++i)
361  {
362  for (j = 0; j < m_expDim; ++j, ++l)
363  {
364  for (k = 0; k < m_coordDim; ++k)
365  {
366  Vmath::Vvtvp(ptsTgt, &deriv[i][k][0], 1,
367  &deriv[j][k][0], 1,
368  &tmp[l][0], 1,
369  &tmp[l][0], 1);
370  }
371  }
372  }
373 
374  Adjoint(tmp, gmat);
375 
376  // Compute g = det(g_{ij}) (= Jacobian squared) and store
377  // temporarily in m_jac.
378  for (i = 0; i < m_expDim; ++i)
379  {
380  Vmath::Vvtvp(ptsTgt, &tmp[i][0], 1, &gmat[i*m_expDim][0], 1,
381  &jac[0], 1, &jac[0], 1);
382  }
383 
384  for (i = 0; i < m_expDim*m_expDim; ++i)
385  {
386  Vmath::Vdiv(ptsTgt, &gmat[i][0], 1, &jac[0], 1, &gmat[i][0], 1);
387  }
388 
389  return gmat;
390  }
391 
392 
393  /**
394  * @param keyTgt Target point distributions.
395  * @returns Derivative factors evaluated at the target
396  * point distributions.
397  */
398  Array<TwoD, NekDouble> GeomFactors::ComputeDerivFactors(
399  const LibUtilities::PointsKeyVector& keyTgt) const
400  {
401  ASSERTL1(keyTgt.size() == m_expDim,
402  "Dimension of target point distribution does not match "
403  "expansion dimension.");
404 
405  int i = 0, j = 0, k = 0, l = 0;
406  int ptsTgt = 1;
407 
408  if (m_type == eDeformed)
409  {
410  // Allocate storage and compute number of points
411  for (i = 0; i < m_expDim; ++i)
412  {
413  ptsTgt *= keyTgt[i].GetNumPoints();
414  }
415  }
416 
417  // Get derivative at geometry points
418  DerivStorage deriv = ComputeDeriv(keyTgt);
419 
420  Array<TwoD, NekDouble> tmp (m_expDim*m_expDim, ptsTgt, 0.0);
421  Array<TwoD, NekDouble> gmat(m_expDim*m_expDim, ptsTgt, 0.0);
422  Array<OneD, NekDouble> jac (ptsTgt, 0.0);
423  Array<TwoD, NekDouble> factors(m_expDim*m_coordDim, ptsTgt, 0.0);
424 
425  // Compute g_{ij} as t_i \cdot t_j and store in tmp
426  for (i = 0, l = 0; i < m_expDim; ++i)
427  {
428  for (j = 0; j < m_expDim; ++j, ++l)
429  {
430  for (k = 0; k < m_coordDim; ++k)
431  {
432  Vmath::Vvtvp(ptsTgt, &deriv[i][k][0], 1,
433  &deriv[j][k][0], 1,
434  &tmp[l][0], 1,
435  &tmp[l][0], 1);
436  }
437  }
438  }
439 
440  Adjoint(tmp, gmat);
441 
442  // Compute g = det(g_{ij}) (= Jacobian squared) and store
443  // temporarily in m_jac.
444  for (i = 0; i < m_expDim; ++i)
445  {
446  Vmath::Vvtvp(ptsTgt, &tmp[i][0], 1, &gmat[i*m_expDim][0], 1,
447  &jac[0], 1, &jac[0], 1);
448  }
449 
450  for (i = 0; i < m_expDim*m_expDim; ++i)
451  {
452  Vmath::Vdiv(ptsTgt, &gmat[i][0], 1, &jac[0], 1, &gmat[i][0], 1);
453  }
454 
455  // Compute the Jacobian = sqrt(g)
456  Vmath::Vsqrt(ptsTgt, &jac[0], 1, &jac[0], 1);
457 
458  // Compute the derivative factors
459  for (k = 0, l = 0; k < m_coordDim; ++k)
460  {
461  for (j = 0; j < m_expDim; ++j, ++l)
462  {
463  for (i = 0; i < m_expDim; ++i)
464  {
465  Vmath::Vvtvp(ptsTgt, &deriv[i][k][0], 1,
466  &gmat[m_expDim*i+j][0], 1,
467  &factors[l][0], 1,
468  &factors[l][0], 1);
469  }
470  }
471  }
472 
473  return factors;
474  }
475 
476 
477  /**
478  * Constructs the Jacobian as per Spencer's book p158 and tests if
479  * negative.
480  */
481  void GeomFactors::CheckIfValid()
482  {
483  // Jacobian test only makes sense when expdim = coorddim
484  // If one-dimensional then element is valid.
485  if (m_coordDim != m_expDim || m_expDim == 1)
486  {
487  m_valid = true;
488  return;
489  }
490 
491  LibUtilities::PointsKeyVector p(m_expDim);
492  int nqtot = 1;
493  for (int i = 0; i < m_expDim; ++i)
494  {
495  p[i] = m_xmap->GetBasis(i)->GetPointsKey();
496  nqtot *= p[i].GetNumPoints();
497  }
498  int pts = (m_type == eRegular || m_type == eMovingRegular)
499  ? 1 : nqtot;
500  Array<OneD, NekDouble> jac(pts, 0.0);
501 
502  DerivStorage deriv = GetDeriv(p);
503 
504  switch (m_expDim)
505  {
506  case 2:
507  {
508  Vmath::Vvtvvtm(pts, &deriv[0][0][0], 1, &deriv[1][1][0], 1,
509  &deriv[1][0][0], 1, &deriv[0][1][0], 1,
510  &jac[0], 1);
511  break;
512  }
513  case 3:
514  {
515  Array<OneD, NekDouble> tmp(pts, 0.0);
516 
517  Vmath::Vvtvvtm(pts, &deriv[1][1][0], 1, &deriv[2][2][0], 1,
518  &deriv[2][1][0], 1, &deriv[1][2][0], 1,
519  &tmp[0], 1);
520  Vmath::Vvtvp (pts, &deriv[0][0][0], 1, &tmp[0], 1,
521  &jac[0], 1, &jac[0], 1);
522 
523  Vmath::Vvtvvtm(pts, &deriv[2][1][0], 1, &deriv[0][2][0], 1,
524  &deriv[0][1][0], 1, &deriv[2][2][0], 1,
525  &tmp[0], 1);
526  Vmath::Vvtvp (pts, &deriv[1][0][0], 1, &tmp[0], 1,
527  &jac[0], 1, &jac[0], 1);
528 
529  Vmath::Vvtvvtm(pts, &deriv[0][1][0], 1, &deriv[1][2][0], 1,
530  &deriv[1][1][0], 1, &deriv[0][2][0], 1,
531  &tmp[0], 1);
532  Vmath::Vvtvp (pts, &deriv[2][0][0], 1, &tmp[0], 1,
533  &jac[0], 1, &jac[0], 1);
534 
535  break;
536  }
537  }
538 
539  if (Vmath::Vmin(pts, &jac[0], 1) < 0)
540  {
541  m_valid = false;
542  }
543  }
544 
545 
546  /**
547  * @param map_points Source data point distribution.
548  * @param src Source data to be interpolated.
549  * @param tpoints Target data point distribution.
550  * @param tgt Target data storage.
551  */
552  void GeomFactors::Interp(
553  const LibUtilities::PointsKeyVector &src_points,
554  const Array<OneD, const NekDouble> &src,
555  const LibUtilities::PointsKeyVector &tgt_points,
556  Array<OneD, NekDouble> &tgt) const
557  {
558  ASSERTL1(src_points.size() == tgt_points.size(),
559  "Dimension of target point distribution does not match "
560  "expansion dimension.");
561 
562  switch (m_expDim)
563  {
564  case 1:
565  LibUtilities::Interp1D(src_points[0], src,
566  tgt_points[0], tgt);
567  break;
568  case 2:
569  LibUtilities::Interp2D(src_points[0], src_points[1], src,
570  tgt_points[0], tgt_points[1], tgt);
571  break;
572  case 3:
573  LibUtilities::Interp3D(src_points[0], src_points[1],
574  src_points[2], src,
575  tgt_points[0], tgt_points[1],
576  tgt_points[2], tgt);
577  break;
578  }
579  }
580 
581 
582  /**
583  * Input and output arrays are of dimension
584  * (m_expDim*m_expDim) x num_points. The first index of the input and
585  * output arrays are ordered row-by-row.
586  * @param src Input data array.
587  * @param tgt Storage for adjoint matrix data.
588  */
589  void GeomFactors::Adjoint(
590  const Array<TwoD, const NekDouble>& src,
591  Array<TwoD, NekDouble>& tgt) const
592  {
593  ASSERTL1(src.num_elements() == tgt.num_elements(),
594  "Source matrix is of different size to destination"
595  "matrix for computing adjoint.");
596 
597  int n = src[0].num_elements();
598  switch (m_expDim)
599  {
600  case 1:
601  Vmath::Fill (n, 1.0, &tgt[0][0], 1);
602  break;
603  case 2:
604  Vmath::Vcopy(n, &src[3][0], 1, &tgt[0][0], 1);
605  Vmath::Smul (n, -1.0, &src[1][0], 1, &tgt[1][0], 1);
606  Vmath::Smul (n, -1.0, &src[2][0], 1, &tgt[2][0], 1);
607  Vmath::Vcopy(n, &src[0][0], 1, &tgt[3][0], 1);
608  break;
609  case 3:
610  {
611  int a, b, c, d, e, i, j;
612 
613  // Compute g^{ij} by computing Cofactors(g_ij)^T
614  for (i = 0; i < m_expDim; ++i)
615  {
616  for (j = 0; j < m_expDim; ++j)
617  {
618  a = ((i+1)%m_expDim) * m_expDim + ((j+1)%m_expDim);
619  b = ((i+1)%m_expDim) * m_expDim + ((j+2)%m_expDim);
620  c = ((i+2)%m_expDim) * m_expDim + ((j+1)%m_expDim);
621  d = ((i+2)%m_expDim) * m_expDim + ((j+2)%m_expDim);
622  e = j*m_expDim + i;
623  Vmath::Vvtvvtm(n, &src[a][0], 1, &src[d][0], 1,
624  &src[b][0], 1, &src[c][0], 1,
625  &tgt[e][0], 1);
626  }
627  }
628  break;
629  }
630  }
631  }
632 
633  }; //end of namespace
634 }; //end of namespace
Generated on Sun Mar 15 2015 00:11:31 for Nektar++ by   doxygen 1.8.1.2