Reference documentation for deal.II version GIT 194dd8bb02 2022-12-03 08:20:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
point.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2021 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 
23 #include <deal.II/base/tensor.h>
24 
26 #include <boost/geometry/core/cs.hpp>
27 #include <boost/geometry/geometries/point.hpp>
29 
30 #include <cmath>
31 
33 
109 template <int dim, typename Number = double>
110 class Point : public Tensor<1, dim, Number>
111 {
112 public:
120  Point();
121 
125  explicit DEAL_II_HOST_DEVICE
127 
136  explicit DEAL_II_HOST_DEVICE
137  Point(const Number x);
138 
149  Point(const Number x, const Number y);
150 
161  Point(const Number x, const Number y, const Number z);
162 
166  template <std::size_t dummy_dim,
167  std::enable_if_t<(dim == dummy_dim) && (dummy_dim != 0), int> = 0>
168  Point(const boost::geometry::model::
169  point<Number, dummy_dim, boost::geometry::cs::cartesian> &boost_pt);
170 
179  unit_vector(const unsigned int i);
180 
186  DEAL_II_HOST_DEVICE Number
187  operator()(const unsigned int index) const;
188 
194  DEAL_II_HOST_DEVICE Number &
195  operator()(const unsigned int index);
196 
202  template <typename OtherNumber>
205 
218 
230 
241 
248  operator-() const;
249 
266  template <typename OtherNumber>
268  dim,
269  typename ProductType<Number,
270  typename EnableIfScalar<OtherNumber>::type>::type>
271  operator*(const OtherNumber) const;
272 
278  template <typename OtherNumber>
280  dim,
281  typename ProductType<Number,
282  typename EnableIfScalar<OtherNumber>::type>::type>
283  operator/(const OtherNumber) const;
284 
290  DEAL_II_HOST_DEVICE Number
292 
306  square() const;
307 
316  distance(const Point<dim, Number> &p) const;
317 
326 
336  template <class Archive>
337  void
338  serialize(Archive &ar, const unsigned int version);
339 };
340 
341 /*--------------------------- Inline functions: Point -----------------------*/
342 
343 #ifndef DOXYGEN
344 
345 // At least clang-3.7 requires us to have a user-defined constructor
346 // and we can't use 'Point<dim,Number>::Point () = default' here.
347 template <int dim, typename Number>
348 inline DEAL_II_HOST_DEVICE
349 Point<dim, Number>::Point() // NOLINT
350 {}
351 
352 
353 
354 template <int dim, typename Number>
355 inline DEAL_II_HOST_DEVICE
357  : Tensor<1, dim, Number>(t)
358 {}
359 
360 
361 
362 template <int dim, typename Number>
363 inline DEAL_II_HOST_DEVICE
364 Point<dim, Number>::Point(const Number x)
365 {
366 # ifndef __CUDA_ARCH__
367  Assert(dim == 1,
368  ExcMessage(
369  "You can only initialize Point<1> objects using the constructor "
370  "that takes only one argument. Point<dim> objects with dim!=1 "
371  "require initialization with the constructor that takes 'dim' "
372  "arguments."));
373 # endif
374 
375  // we can only get here if we pass the assertion. use the switch anyway so
376  // as to avoid compiler warnings about uninitialized elements or writing
377  // beyond the end of the 'values' array
378  switch (dim)
379  {
380  case 1:
381  this->values[0] = x;
382  break;
383 
384  default:;
385  }
386 }
387 
388 
389 
390 template <int dim, typename Number>
391 inline DEAL_II_HOST_DEVICE
392 Point<dim, Number>::Point(const Number x, const Number y)
393 {
394 # ifndef __CUDA_ARCH__
395  Assert(dim == 2,
396  ExcMessage(
397  "You can only initialize Point<2> objects using the constructor "
398  "that takes two arguments. Point<dim> objects with dim!=2 "
399  "require initialization with the constructor that takes 'dim' "
400  "arguments."));
401 # endif
402 
403  // we can only get here if we pass the assertion. use the indirection anyway
404  // so as to avoid compiler warnings about uninitialized elements or writing
405  // beyond the end of the 'values' array
406  constexpr unsigned int y_index = (dim < 2) ? 0 : 1;
407  this->values[0] = x;
408  this->values[y_index] = y;
409 }
410 
411 
412 
413 template <int dim, typename Number>
414 inline DEAL_II_HOST_DEVICE
415 Point<dim, Number>::Point(const Number x, const Number y, const Number z)
416 {
417 # ifndef __CUDA_ARCH__
418  Assert(dim == 3,
419  ExcMessage(
420  "You can only initialize Point<3> objects using the constructor "
421  "that takes three arguments. Point<dim> objects with dim!=3 "
422  "require initialization with the constructor that takes 'dim' "
423  "arguments."));
424 # endif
425 
426  // we can only get here if we pass the assertion. use the indirection anyway
427  // so as to avoid compiler warnings about uninitialized elements or writing
428  // beyond the end of the 'values' array
429  constexpr unsigned int y_index = (dim < 2) ? 0 : 1;
430  constexpr unsigned int z_index = (dim < 3) ? 0 : 2;
431  this->values[0] = x;
432  this->values[y_index] = y;
433  this->values[z_index] = z;
434 }
435 
436 
437 
438 template <int dim, typename Number>
439 template <std::size_t dummy_dim,
440  std::enable_if_t<(dim == dummy_dim) && (dummy_dim != 0), int>>
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>
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_HOST_DEVICE Number
471 Point<dim, Number>::operator()(const unsigned int index) const
472 {
473 # ifndef __CUDA_ARCH__
474  AssertIndexRange(static_cast<int>(index), dim);
475 # endif
476  return this->values[index];
477 }
478 
479 
480 
481 template <int dim, typename Number>
482 inline DEAL_II_HOST_DEVICE Number &
483 Point<dim, Number>::operator()(const unsigned int index)
484 {
485 # ifndef __CUDA_ARCH__
486  AssertIndexRange(static_cast<int>(index), dim);
487 # endif
488  return this->values[index];
489 }
490 
491 
492 
493 template <int dim, typename Number>
494 template <typename OtherNumber>
497 {
499  return *this;
500 }
501 
502 
503 
504 template <int dim, typename Number>
507 {
508  Point<dim, Number> tmp = *this;
509  tmp += p;
510  return tmp;
511 }
512 
513 
514 
515 template <int dim, typename Number>
518 {
519  return (Tensor<1, dim, Number>(*this) -= p);
520 }
521 
522 
523 
524 template <int dim, typename Number>
527 {
528  Point<dim, Number> tmp = *this;
529  tmp -= p;
530  return tmp;
531 }
532 
533 
534 
535 template <int dim, typename 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_HOST_DEVICE
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_HOST_DEVICE
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_HOST_DEVICE Number
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>
595 {
596  return this->norm_square();
597 }
598 
599 
600 
601 template <int dim, typename Number>
604 {
605  return std::sqrt(distance_square(p));
606 }
607 
608 
609 
610 template <int dim, typename Number>
613 {
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 
649 template <int dim, typename Number, typename OtherNumber>
650 inline DEAL_II_HOST_DEVICE
651  Point<dim,
652  typename ProductType<Number,
653  typename EnableIfScalar<OtherNumber>::type>::type>
654  operator*(const OtherNumber factor, const Point<dim, Number> &p)
655 {
656  return p * factor;
657 }
658 
659 
660 
666 template <int dim, typename Number>
667 inline std::ostream &
668 operator<<(std::ostream &out, const Point<dim, Number> &p)
669 {
670  for (unsigned int i = 0; i < dim - 1; ++i)
671  out << p[i] << ' ';
672  out << p[dim - 1];
673 
674  return out;
675 }
676 
677 
678 
683 template <int dim, typename Number>
684 inline std::istream &
685 operator>>(std::istream &in, Point<dim, Number> &p)
686 {
687  for (unsigned int i = 0; i < dim; ++i)
688  in >> p[i];
689 
690  return in;
691 }
692 
693 
694 #ifndef DOXYGEN
695 
701 template <typename Number>
702 inline std::ostream &
703 operator<<(std::ostream &out, const Point<1, Number> &p)
704 {
705  out << p[0];
706 
707  return out;
708 }
709 
710 #endif // DOXYGEN
712 
713 #endif
OutputOperator< VectorType > & operator<<(OutputOperator< VectorType > &out, unsigned int step)
Definition: operator.h:165
Definition: point.h:111
numbers::NumberTraits< Number >::real_type distance_square(const Point< dim, Number > &p) const
Point< dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const OtherNumber) const
Tensor< 1, dim, Number > operator-(const Point< dim, Number > &) const
Point< dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const OtherNumber) const
Number operator*(const Tensor< 1, dim, Number > &p) const
Number & operator()(const unsigned int index)
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
Point(const boost::geometry::model::point< Number, dummy_dim, boost::geometry::cs::cartesian > &boost_pt)
Point(const Tensor< 1, dim, Number > &)
Point(const Number x, const Number y, const Number z)
Point< dim, Number > & operator=(const Tensor< 1, dim, OtherNumber > &p)
Point< dim, Number > operator-(const Tensor< 1, dim, Number > &) const
Point(const Number x)
void serialize(Archive &ar, const unsigned int version)
numbers::NumberTraits< Number >::real_type square() const
std::istream & operator>>(std::istream &in, Point< dim, Number > &p)
Definition: point.h:685
Point< dim, Number > operator-() const
Point< dim, Number > operator+(const Tensor< 1, dim, Number > &) const
Point(const Number x, const Number y)
static Point< dim, Number > unit_vector(const unsigned int i)
Number operator()(const unsigned int index) const
Point< dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const OtherNumber factor, const Point< dim, Number > &p)
Definition: point.h:654
Definition: tensor.h:503
constexpr Tensor & operator=(const Tensor< rank_, dim, OtherNumber > &rhs)
constexpr numbers::NumberTraits< double >::real_type norm_square() const
Tensor< rank_ - 1, dim, double > values[(dim !=0) ? dim :1]
Definition: tensor.h:864
Tensor< rank, dim, double > sum(const Tensor< rank, dim, double > &local, const MPI_Comm &mpi_communicator)
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:102
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:512
#define Assert(cond, exc)
Definition: exceptions.h:1501
static ::ExceptionBase & ExcNotImplemented()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1760
static ::ExceptionBase & ExcMessage(std::string arg1)
#define DEAL_II_HOST_DEVICE
Definition: numbers.h:35
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T & value(const T &t)
Definition: numbers.h:733
static constexpr std::enable_if_t< std::is_same< Dummy, number >::value &&is_cuda_compatible< Dummy >::value, real_type > abs_square(const number &x)