Reference documentation for deal.II version Git 71593a1301 2019-07-19 12:10:10 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
point.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_point_h
17 #define dealii_point_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/exceptions.h>
23 #include <deal.II/base/tensor.h>
24 
25 DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
26 #include <boost/geometry.hpp>
27 DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
28 
29 #include <cmath>
30 
31 DEAL_II_NAMESPACE_OPEN
32 
109 template <int dim, typename Number = double>
110 class Point : public Tensor<1, dim, Number>
111 {
112 public:
119  DEAL_II_CUDA_HOST_DEV
120  Point();
121 
125  explicit DEAL_II_CUDA_HOST_DEV
126  Point(const Tensor<1, dim, Number> &);
127 
136  explicit DEAL_II_CUDA_HOST_DEV
137  Point(const Number x);
138 
148  DEAL_II_CUDA_HOST_DEV
149  Point(const Number x, const Number y);
150 
160  DEAL_II_CUDA_HOST_DEV
161  Point(const Number x, const Number y, const Number z);
162 
166  template <std::size_t dummy_dim,
167  typename std::enable_if<(dim == dummy_dim) && (dummy_dim != 0),
168  int>::type = 0>
169  Point(const boost::geometry::model::
170  point<Number, dummy_dim, boost::geometry::cs::cartesian> &boost_pt);
171 
179  static DEAL_II_CUDA_HOST_DEV Point<dim, Number>
180  unit_vector(const unsigned int i);
181 
187  DEAL_II_CUDA_HOST_DEV Number
188  operator()(const unsigned int index) const;
189 
195  DEAL_II_CUDA_HOST_DEV Number &
196  operator()(const unsigned int index);
197 
208  DEAL_II_CUDA_HOST_DEV Point<dim, Number>
209  operator+(const Tensor<1, dim, Number> &) const;
210 
220  DEAL_II_CUDA_HOST_DEV Tensor<1, dim, Number>
221  operator-(const Point<dim, Number> &) const;
222 
231  DEAL_II_CUDA_HOST_DEV Point<dim, Number>
232  operator-(const Tensor<1, dim, Number> &) const;
233 
239  DEAL_II_CUDA_HOST_DEV Point<dim, Number>
240  operator-() const;
241 
258  template <typename OtherNumber>
259  DEAL_II_CUDA_HOST_DEV Point<
260  dim,
261  typename ProductType<Number,
262  typename EnableIfScalar<OtherNumber>::type>::type>
263  operator*(const OtherNumber) const;
264 
270  template <typename OtherNumber>
271  DEAL_II_CUDA_HOST_DEV Point<
272  dim,
273  typename ProductType<Number,
274  typename EnableIfScalar<OtherNumber>::type>::type>
275  operator/(const OtherNumber) const;
276 
282  DEAL_II_CUDA_HOST_DEV Number operator*(const Tensor<1, dim, Number> &p) const;
283 
296  DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
297  square() const;
298 
306  DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
307  distance(const Point<dim, Number> &p) const;
308 
315  DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
316  distance_square(const Point<dim, Number> &p) const;
317 
326  template <class Archive>
327  void
328  serialize(Archive &ar, const unsigned int version);
329 };
330 
331 /*--------------------------- Inline functions: Point -----------------------*/
332 
333 #ifndef DOXYGEN
334 
335 // At least clang-3.7 requires us to have a user-defined constructor
336 // and we can't use 'Point<dim,Number>::Point () = default' here.
337 template <int dim, typename Number>
338 inline DEAL_II_CUDA_HOST_DEV
339 Point<dim, Number>::Point() // NOLINT
340 {}
341 
342 
343 
344 template <int dim, typename Number>
345 inline DEAL_II_CUDA_HOST_DEV
348 {}
349 
350 
351 
352 template <int dim, typename Number>
353 inline DEAL_II_CUDA_HOST_DEV
354 Point<dim, Number>::Point(const Number x)
355 {
356 # ifndef __CUDA_ARCH__
357  Assert(dim == 1,
358  ExcMessage(
359  "You can only initialize Point<1> objects using the constructor "
360  "that takes only one argument. Point<dim> objects with dim!=1 "
361  "require initialization with the constructor that takes 'dim' "
362  "arguments."));
363 # endif
364 
365  // we can only get here if we pass the assertion. use the switch anyway so
366  // as to avoid compiler warnings about uninitialized elements or writing
367  // beyond the end of the 'values' array
368  switch (dim)
369  {
370  case 1:
371  this->values[0] = x;
372  break;
373 
374  default:;
375  }
376 }
377 
378 
379 
380 template <int dim, typename Number>
381 inline DEAL_II_CUDA_HOST_DEV
382 Point<dim, Number>::Point(const Number x, const Number y)
383 {
384 # ifndef __CUDA_ARCH__
385  Assert(dim == 2,
386  ExcMessage(
387  "You can only initialize Point<2> objects using the constructor "
388  "that takes two arguments. Point<dim> objects with dim!=2 "
389  "require initialization with the constructor that takes 'dim' "
390  "arguments."));
391 # endif
392 
393  // we can only get here if we pass the assertion. use the indirection anyway
394  // so as to avoid compiler warnings about uninitialized elements or writing
395  // beyond the end of the 'values' array
396  constexpr unsigned int y_index = (dim < 2) ? 0 : 1;
397  this->values[0] = x;
398  this->values[y_index] = y;
399 }
400 
401 
402 
403 template <int dim, typename Number>
404 inline DEAL_II_CUDA_HOST_DEV
405 Point<dim, Number>::Point(const Number x, const Number y, const Number z)
406 {
407 # ifndef __CUDA_ARCH__
408  Assert(dim == 3,
409  ExcMessage(
410  "You can only initialize Point<3> objects using the constructor "
411  "that takes three arguments. Point<dim> objects with dim!=3 "
412  "require initialization with the constructor that takes 'dim' "
413  "arguments."));
414 # endif
415 
416  // we can only get here if we pass the assertion. use the indirection anyway
417  // so as to avoid compiler warnings about uninitialized elements or writing
418  // beyond the end of the 'values' array
419  constexpr unsigned int y_index = (dim < 2) ? 0 : 1;
420  constexpr unsigned int z_index = (dim < 3) ? 0 : 2;
421  this->values[0] = x;
422  this->values[y_index] = y;
423  this->values[z_index] = z;
424 }
425 
426 
427 
428 template <int dim, typename Number>
429 template <
430  std::size_t dummy_dim,
431  typename std::enable_if<(dim == dummy_dim) && (dummy_dim != 0), int>::type>
433  const boost::geometry::model::
434  point<Number, dummy_dim, boost::geometry::cs::cartesian> &boost_pt)
435 {
436  Assert(dim <= 3, ExcNotImplemented());
437  this->values[0] = boost::geometry::get<0>(boost_pt);
438  constexpr unsigned int y_index = (dim < 2) ? 0 : 1;
439  constexpr unsigned int z_index = (dim < 3) ? 0 : 2;
440 
441  if (dim >= 2)
442  this->values[y_index] = boost::geometry::get<y_index>(boost_pt);
443 
444  if (dim >= 3)
445  this->values[z_index] = boost::geometry::get<z_index>(boost_pt);
446 }
447 
448 
449 
450 template <int dim, typename Number>
451 inline DEAL_II_CUDA_HOST_DEV Point<dim, Number>
452  Point<dim, Number>::unit_vector(unsigned int i)
453 {
455  p[i] = 1.;
456  return p;
457 }
458 
459 
460 template <int dim, typename Number>
461 inline DEAL_II_CUDA_HOST_DEV Number
462 Point<dim, Number>::operator()(const unsigned int index) const
463 {
464 # ifndef __CUDA_ARCH__
465  AssertIndexRange(index, dim);
466 # endif
467  return this->values[index];
468 }
469 
470 
471 
472 template <int dim, typename Number>
473 inline DEAL_II_CUDA_HOST_DEV Number &
474 Point<dim, Number>::operator()(const unsigned int index)
475 {
476 # ifndef __CUDA_ARCH__
477  AssertIndexRange(index, dim);
478 # endif
479  return this->values[index];
480 }
481 
482 
483 
484 template <int dim, typename Number>
485 inline DEAL_II_CUDA_HOST_DEV Point<dim, Number>
487 {
488  Point<dim, Number> tmp = *this;
489  tmp += p;
490  return tmp;
491 }
492 
493 
494 
495 template <int dim, typename Number>
496 inline DEAL_II_CUDA_HOST_DEV Tensor<1, dim, Number>
498 {
499  return (Tensor<1, dim, Number>(*this) -= p);
500 }
501 
502 
503 
504 template <int dim, typename Number>
505 inline DEAL_II_CUDA_HOST_DEV Point<dim, Number>
507 {
508  Point<dim, Number> tmp = *this;
509  tmp -= p;
510  return tmp;
511 }
512 
513 
514 
515 template <int dim, typename Number>
516 inline DEAL_II_CUDA_HOST_DEV Point<dim, Number>
518 {
519  Point<dim, Number> result;
520  for (unsigned int i = 0; i < dim; ++i)
521  result.values[i] = -this->values[i];
522  return result;
523 }
524 
525 
526 
527 template <int dim, typename Number>
528 template <typename OtherNumber>
529 inline DEAL_II_CUDA_HOST_DEV
530  Point<dim,
531  typename ProductType<Number,
532  typename EnableIfScalar<OtherNumber>::type>::type>
533  Point<dim, Number>::operator*(const OtherNumber factor) const
534 {
536  for (unsigned int i = 0; i < dim; ++i)
537  tmp[i] = this->operator[](i) * factor;
538  return tmp;
539 }
540 
541 
542 
543 template <int dim, typename Number>
544 template <typename OtherNumber>
545 inline DEAL_II_CUDA_HOST_DEV
546  Point<dim,
547  typename ProductType<Number,
548  typename EnableIfScalar<OtherNumber>::type>::type>
549  Point<dim, Number>::operator/(const OtherNumber factor) const
550 {
552  for (unsigned int i = 0; i < dim; ++i)
553  tmp[i] = this->operator[](i) / factor;
554  return tmp;
555 }
556 
557 
558 
559 template <int dim, typename Number>
560 inline DEAL_II_CUDA_HOST_DEV Number Point<dim, Number>::
561  operator*(const Tensor<1, dim, Number> &p) const
562 {
563  Number res = Number();
564  for (unsigned int i = 0; i < dim; ++i)
565  res += this->operator[](i) * p[i];
566  return res;
567 }
568 
569 
570 template <int dim, typename Number>
571 inline DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
573 {
574  return this->norm_square();
575 }
576 
577 
578 
579 template <int dim, typename Number>
580 inline DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
582 {
583  return std::sqrt(distance_square(p));
584 }
585 
586 
587 
588 template <int dim, typename Number>
589 inline DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
591 {
593  for (unsigned int i = 0; i < dim; ++i)
594  {
595  const Number diff = static_cast<Number>(this->values[i]) - p(i);
597  }
598 
599  return sum;
600 }
601 
602 
603 
604 template <int dim, typename Number>
605 template <class Archive>
606 inline void
607 Point<dim, Number>::serialize(Archive &ar, const unsigned int)
608 {
609  // forward to serialization
610  // function in the base class
611  ar &static_cast<Tensor<1, dim, Number> &>(*this);
612 }
613 
614 #endif // DOXYGEN
615 
616 
617 /*--------------------------- Global functions: Point -----------------------*/
618 
619 
628 template <int dim, typename Number, typename OtherNumber>
629 inline DEAL_II_CUDA_HOST_DEV
630  Point<dim,
631  typename ProductType<Number,
632  typename EnableIfScalar<OtherNumber>::type>::type>
633  operator*(const OtherNumber factor, const Point<dim, Number> &p)
634 {
635  return p * factor;
636 }
637 
638 
639 
645 template <int dim, typename Number>
646 inline std::ostream &
647 operator<<(std::ostream &out, const Point<dim, Number> &p)
648 {
649  for (unsigned int i = 0; i < dim - 1; ++i)
650  out << p[i] << ' ';
651  out << p[dim - 1];
652 
653  return out;
654 }
655 
656 
657 
662 template <int dim, typename Number>
663 inline std::istream &
664 operator>>(std::istream &in, Point<dim, Number> &p)
665 {
666  for (unsigned int i = 0; i < dim; ++i)
667  in >> p[i];
668 
669  return in;
670 }
671 
672 
673 #ifndef DOXYGEN
674 
680 template <typename Number>
681 inline std::ostream &
682 operator<<(std::ostream &out, const Point<1, Number> &p)
683 {
684  out << p[0];
685 
686  return out;
687 }
688 
689 #endif // DOXYGEN
690 DEAL_II_NAMESPACE_CLOSE
691 
692 #endif
Tensor< rank, dim, Number > sum(const Tensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
Point< dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const OtherNumber factor, const Point< dim, Number > &p)
Definition: point.h:633
Point< dim, Number > operator+(const Tensor< 1, dim, Number > &) const
static constexpr std::enable_if< std::is_same< Dummy, number >::value &&is_cuda_compatible< Dummy >::value, real_type >::type abs_square(const number &x)
Definition: numbers.h:560
#define AssertIndexRange(index, range)
Definition: exceptions.h:1637
static Point< dim, Number > unit_vector(const unsigned int i)
Tensor< rank_ - 1, dim, Number > values[(dim !=0) ? dim :1]
Definition: tensor.h:713
Number operator()(const unsigned int index) const
Definition: point.h:110
static ::ExceptionBase & ExcMessage(std::string arg1)
Point< dim, Number > operator-() const
#define Assert(cond, exc)
Definition: exceptions.h:1407
Point< dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const OtherNumber) const
void serialize(Archive &ar, const unsigned int version)
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
DEAL_II_CONSTEXPR numbers::NumberTraits< Number >::real_type norm_square() const
Definition: tensor.h:1395
Definition: mpi.h:91
numbers::NumberTraits< Number >::real_type square() const
Point< dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const OtherNumber) const
static ::ExceptionBase & ExcNotImplemented()
numbers::NumberTraits< Number >::real_type distance_square(const Point< dim, Number > &p) const
std::istream & operator>>(std::istream &in, Point< dim, Number > &p)
Definition: point.h:664