Reference documentation for deal.II version Git 0491297983 2019-09-23 09:31:59 +0200
\(\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 
203  template <typename OtherNumber>
206 
217  DEAL_II_CUDA_HOST_DEV Point<dim, Number>
218  operator+(const Tensor<1, dim, Number> &) const;
219 
229  DEAL_II_CUDA_HOST_DEV Tensor<1, dim, Number>
230  operator-(const Point<dim, Number> &) const;
231 
240  DEAL_II_CUDA_HOST_DEV Point<dim, Number>
241  operator-(const Tensor<1, dim, Number> &) const;
242 
248  DEAL_II_CUDA_HOST_DEV Point<dim, Number>
249  operator-() const;
250 
267  template <typename OtherNumber>
268  DEAL_II_CUDA_HOST_DEV Point<
269  dim,
270  typename ProductType<Number,
271  typename EnableIfScalar<OtherNumber>::type>::type>
272  operator*(const OtherNumber) const;
273 
279  template <typename OtherNumber>
280  DEAL_II_CUDA_HOST_DEV Point<
281  dim,
282  typename ProductType<Number,
283  typename EnableIfScalar<OtherNumber>::type>::type>
284  operator/(const OtherNumber) const;
285 
291  DEAL_II_CUDA_HOST_DEV Number operator*(const Tensor<1, dim, Number> &p) const;
292 
305  DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
306  square() const;
307 
315  DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
316  distance(const Point<dim, Number> &p) const;
317 
324  DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
325  distance_square(const Point<dim, Number> &p) const;
326 
335  template <class Archive>
336  void
337  serialize(Archive &ar, const unsigned int version);
338 };
339 
340 /*--------------------------- Inline functions: Point -----------------------*/
341 
342 #ifndef DOXYGEN
343 
344 // At least clang-3.7 requires us to have a user-defined constructor
345 // and we can't use 'Point<dim,Number>::Point () = default' here.
346 template <int dim, typename Number>
347 inline DEAL_II_CUDA_HOST_DEV
348 Point<dim, Number>::Point() // NOLINT
349 {}
350 
351 
352 
353 template <int dim, typename Number>
354 inline DEAL_II_CUDA_HOST_DEV
357 {}
358 
359 
360 
361 template <int dim, typename Number>
362 inline DEAL_II_CUDA_HOST_DEV
363 Point<dim, Number>::Point(const Number x)
364 {
365 # ifndef __CUDA_ARCH__
366  Assert(dim == 1,
367  ExcMessage(
368  "You can only initialize Point<1> objects using the constructor "
369  "that takes only one argument. Point<dim> objects with dim!=1 "
370  "require initialization with the constructor that takes 'dim' "
371  "arguments."));
372 # endif
373 
374  // we can only get here if we pass the assertion. use the switch anyway so
375  // as to avoid compiler warnings about uninitialized elements or writing
376  // beyond the end of the 'values' array
377  switch (dim)
378  {
379  case 1:
380  this->values[0] = x;
381  break;
382 
383  default:;
384  }
385 }
386 
387 
388 
389 template <int dim, typename Number>
390 inline DEAL_II_CUDA_HOST_DEV
391 Point<dim, Number>::Point(const Number x, const Number y)
392 {
393 # ifndef __CUDA_ARCH__
394  Assert(dim == 2,
395  ExcMessage(
396  "You can only initialize Point<2> objects using the constructor "
397  "that takes two arguments. Point<dim> objects with dim!=2 "
398  "require initialization with the constructor that takes 'dim' "
399  "arguments."));
400 # endif
401 
402  // we can only get here if we pass the assertion. use the indirection anyway
403  // so as to avoid compiler warnings about uninitialized elements or writing
404  // beyond the end of the 'values' array
405  constexpr unsigned int y_index = (dim < 2) ? 0 : 1;
406  this->values[0] = x;
407  this->values[y_index] = y;
408 }
409 
410 
411 
412 template <int dim, typename Number>
413 inline DEAL_II_CUDA_HOST_DEV
414 Point<dim, Number>::Point(const Number x, const Number y, const Number z)
415 {
416 # ifndef __CUDA_ARCH__
417  Assert(dim == 3,
418  ExcMessage(
419  "You can only initialize Point<3> objects using the constructor "
420  "that takes three arguments. Point<dim> objects with dim!=3 "
421  "require initialization with the constructor that takes 'dim' "
422  "arguments."));
423 # endif
424 
425  // we can only get here if we pass the assertion. use the indirection anyway
426  // so as to avoid compiler warnings about uninitialized elements or writing
427  // beyond the end of the 'values' array
428  constexpr unsigned int y_index = (dim < 2) ? 0 : 1;
429  constexpr unsigned int z_index = (dim < 3) ? 0 : 2;
430  this->values[0] = x;
431  this->values[y_index] = y;
432  this->values[z_index] = z;
433 }
434 
435 
436 
437 template <int dim, typename Number>
438 template <
439  std::size_t dummy_dim,
440  typename std::enable_if<(dim == dummy_dim) && (dummy_dim != 0), int>::type>
442  const boost::geometry::model::
443  point<Number, dummy_dim, boost::geometry::cs::cartesian> &boost_pt)
444 {
445  Assert(dim <= 3, ExcNotImplemented());
446  this->values[0] = boost::geometry::get<0>(boost_pt);
447  constexpr unsigned int y_index = (dim < 2) ? 0 : 1;
448  constexpr unsigned int z_index = (dim < 3) ? 0 : 2;
449 
450  if (dim >= 2)
451  this->values[y_index] = boost::geometry::get<y_index>(boost_pt);
452 
453  if (dim >= 3)
454  this->values[z_index] = boost::geometry::get<z_index>(boost_pt);
455 }
456 
457 
458 
459 template <int dim, typename Number>
460 inline DEAL_II_CUDA_HOST_DEV Point<dim, Number>
461  Point<dim, Number>::unit_vector(unsigned int i)
462 {
464  p[i] = 1.;
465  return p;
466 }
467 
468 
469 template <int dim, typename Number>
470 inline DEAL_II_CUDA_HOST_DEV Number
471 Point<dim, Number>::operator()(const unsigned int index) const
472 {
473 # ifndef __CUDA_ARCH__
474  AssertIndexRange(index, dim);
475 # endif
476  return this->values[index];
477 }
478 
479 
480 
481 template <int dim, typename Number>
482 inline DEAL_II_CUDA_HOST_DEV Number &
483 Point<dim, Number>::operator()(const unsigned int index)
484 {
485 # ifndef __CUDA_ARCH__
486  AssertIndexRange(index, dim);
487 # endif
488  return this->values[index];
489 }
490 
491 
492 
493 template <int dim, typename Number>
494 template <typename OtherNumber>
495 inline DEAL_II_ALWAYS_INLINE Point<dim, Number> &
497 {
499  return *this;
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 Tensor<1, dim, Number>
518 {
519  return (Tensor<1, dim, Number>(*this) -= p);
520 }
521 
522 
523 
524 template <int dim, typename Number>
525 inline DEAL_II_CUDA_HOST_DEV Point<dim, Number>
527 {
528  Point<dim, Number> tmp = *this;
529  tmp -= p;
530  return tmp;
531 }
532 
533 
534 
535 template <int dim, typename Number>
536 inline DEAL_II_CUDA_HOST_DEV Point<dim, Number>
538 {
539  Point<dim, Number> result;
540  for (unsigned int i = 0; i < dim; ++i)
541  result.values[i] = -this->values[i];
542  return result;
543 }
544 
545 
546 
547 template <int dim, typename Number>
548 template <typename OtherNumber>
549 inline DEAL_II_CUDA_HOST_DEV
550  Point<dim,
551  typename ProductType<Number,
552  typename EnableIfScalar<OtherNumber>::type>::type>
553  Point<dim, Number>::operator*(const OtherNumber factor) const
554 {
556  for (unsigned int i = 0; i < dim; ++i)
557  tmp[i] = this->operator[](i) * factor;
558  return tmp;
559 }
560 
561 
562 
563 template <int dim, typename Number>
564 template <typename OtherNumber>
565 inline DEAL_II_CUDA_HOST_DEV
566  Point<dim,
567  typename ProductType<Number,
568  typename EnableIfScalar<OtherNumber>::type>::type>
569  Point<dim, Number>::operator/(const OtherNumber factor) const
570 {
571  const Tensor<1, dim, Number> &base_object = *this;
572  return Point<
573  dim,
574  typename ProductType<Number,
575  typename EnableIfScalar<OtherNumber>::type>::type>(
576  ::operator/(base_object, factor));
577 }
578 
579 
580 
581 template <int dim, typename Number>
582 inline DEAL_II_CUDA_HOST_DEV Number Point<dim, Number>::
583  operator*(const Tensor<1, dim, Number> &p) const
584 {
585  Number res = Number();
586  for (unsigned int i = 0; i < dim; ++i)
587  res += this->operator[](i) * p[i];
588  return res;
589 }
590 
591 
592 template <int dim, typename Number>
593 inline DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
595 {
596  return this->norm_square();
597 }
598 
599 
600 
601 template <int dim, typename Number>
602 inline DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
604 {
605  return std::sqrt(distance_square(p));
606 }
607 
608 
609 
610 template <int dim, typename Number>
611 inline DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
613 {
614  Number sum = internal::NumberType<Number>::value(0.0);
615  for (unsigned int i = 0; i < dim; ++i)
616  {
617  const Number diff = static_cast<Number>(this->values[i]) - p(i);
619  }
620 
621  return sum;
622 }
623 
624 
625 
626 template <int dim, typename Number>
627 template <class Archive>
628 inline void
629 Point<dim, Number>::serialize(Archive &ar, const unsigned int)
630 {
631  // forward to serialization
632  // function in the base class
633  ar &static_cast<Tensor<1, dim, Number> &>(*this);
634 }
635 
636 #endif // DOXYGEN
637 
638 
639 /*--------------------------- Global functions: Point -----------------------*/
640 
641 
650 template <int dim, typename Number, typename OtherNumber>
651 inline DEAL_II_CUDA_HOST_DEV
652  Point<dim,
653  typename ProductType<Number,
654  typename EnableIfScalar<OtherNumber>::type>::type>
655  operator*(const OtherNumber factor, const Point<dim, Number> &p)
656 {
657  return p * factor;
658 }
659 
660 
661 
667 template <int dim, typename Number>
668 inline std::ostream &
669 operator<<(std::ostream &out, const Point<dim, Number> &p)
670 {
671  for (unsigned int i = 0; i < dim - 1; ++i)
672  out << p[i] << ' ';
673  out << p[dim - 1];
674 
675  return out;
676 }
677 
678 
679 
684 template <int dim, typename Number>
685 inline std::istream &
686 operator>>(std::istream &in, Point<dim, Number> &p)
687 {
688  for (unsigned int i = 0; i < dim; ++i)
689  in >> p[i];
690 
691  return in;
692 }
693 
694 
695 #ifndef DOXYGEN
696 
702 template <typename Number>
703 inline std::ostream &
704 operator<<(std::ostream &out, const Point<1, Number> &p)
705 {
706  out << p[0];
707 
708  return out;
709 }
710 
711 #endif // DOXYGEN
712 DEAL_II_NAMESPACE_CLOSE
713 
714 #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:655
Point< dim, Number > operator+(const Tensor< 1, dim, Number > &) const
Point< dim, Number > & operator=(const Tensor< 1, dim, OtherNumber > &p)
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:605
#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:718
DEAL_II_CONSTEXPR Tensor & operator=(const Tensor< rank_, dim, OtherNumber > &rhs)
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:1425
Definition: tensor.h:422
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:686