Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-3041-g9c4075ddf4 2025-04-07 16:30:00+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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
point.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1998 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_point_h
16#define dealii_point_h
17
18
19#include <deal.II/base/config.h>
20
24#include <deal.II/base/tensor.h>
25
26#include <boost/geometry/core/cs.hpp>
27#include <boost/geometry/geometries/point.hpp>
28
29#include <cmath>
30
32
110template <int dim, typename Number = double>
112class Point : public Tensor<1, dim, Number>
113{
114public:
121 constexpr DEAL_II_HOST_DEVICE
123
127 constexpr explicit DEAL_II_HOST_DEVICE
129
138 constexpr explicit DEAL_II_HOST_DEVICE
139 Point(const Number x);
140
150 constexpr DEAL_II_HOST_DEVICE
151 Point(const Number x, const Number y);
152
162 constexpr DEAL_II_HOST_DEVICE
163 Point(const Number x, const Number y, const Number z);
164
168 template <std::size_t dummy_dim,
169 std::enable_if_t<(dim == dummy_dim) && (dummy_dim != 0), int> = 0>
170 constexpr Point(
171 const boost::geometry::model::
172 point<Number, dummy_dim, boost::geometry::cs::cartesian> &boost_pt);
173
182 unit_vector(const unsigned int i);
183
189 constexpr DEAL_II_HOST_DEVICE Number
190 operator()(const unsigned int index) const;
191
197 constexpr DEAL_II_HOST_DEVICE Number &
198 operator()(const unsigned int index);
199
205 template <typename OtherNumber>
206 constexpr Point<dim, Number> &
208
221
233
244
251 operator-() const;
252
269 template <typename OtherNumber>
270 constexpr DEAL_II_HOST_DEVICE Point<
271 dim,
272 typename ProductType<Number,
273 typename EnableIfScalar<OtherNumber>::type>::type>
274 operator*(const OtherNumber) const;
275
281 template <typename OtherNumber>
282 constexpr DEAL_II_HOST_DEVICE Point<
283 dim,
284 typename ProductType<Number,
285 typename EnableIfScalar<OtherNumber>::type>::type>
286 operator/(const OtherNumber) const;
287
293 constexpr DEAL_II_HOST_DEVICE Number
295
308 constexpr DEAL_II_HOST_DEVICE
310 square() const;
311
321
328 constexpr DEAL_II_HOST_DEVICE
331
341 template <class Archive>
342 void
343 serialize(Archive &ar, const unsigned int version);
344};
345
346/*--------------------------- Inline functions: Point -----------------------*/
347
348#ifndef DOXYGEN
349
350// At least clang-3.7 requires us to have a user-defined constructor
351// and we can't use 'Point<dim,Number>::Point () = default' here.
352template <int dim, typename Number>
354constexpr DEAL_II_HOST_DEVICE Point<dim, Number>::Point() // NOLINT
355{}
356
357
358
359template <int dim, typename Number>
361constexpr DEAL_II_HOST_DEVICE Point<dim, Number>::Point(
362 const Tensor<1, dim, Number> &t)
363 : Tensor<1, dim, Number>(t)
364{}
365
366
367
368template <int dim, typename Number>
370constexpr DEAL_II_HOST_DEVICE Point<dim, Number>::Point(const Number x)
371{
372 Assert(dim == 1,
374 "You can only initialize Point<1> objects using the constructor "
375 "that takes only one argument. Point<dim> objects with dim!=1 "
376 "require initialization with the constructor that takes 'dim' "
377 "arguments."));
378
379 // we can only get here if we pass the assertion. use the switch anyway so
380 // as to avoid compiler warnings about uninitialized elements or writing
381 // beyond the end of the 'values' array
382 switch (dim)
383 {
384 case 1:
385 this->values[0] = x;
386 break;
387
388 default:;
389 }
390}
391
392
393
394template <int dim, typename Number>
396constexpr DEAL_II_HOST_DEVICE Point<dim, Number>::Point(const Number x,
397 const Number y)
398{
399 Assert(dim == 2,
401 "You can only initialize Point<2> objects using the constructor "
402 "that takes two arguments. Point<dim> objects with dim!=2 "
403 "require initialization with the constructor that takes 'dim' "
404 "arguments."));
405
406 // we can only get here if we pass the assertion. use the indirection anyway
407 // so as to avoid compiler warnings about uninitialized elements or writing
408 // beyond the end of the 'values' array
409 constexpr unsigned int y_index = (dim < 2) ? 0 : 1;
410 this->values[0] = x;
411 this->values[y_index] = y;
412}
413
414
415
416template <int dim, typename Number>
418constexpr DEAL_II_HOST_DEVICE Point<dim, Number>::Point(const Number x,
419 const Number y,
420 const Number z)
421{
422 Assert(dim == 3,
424 "You can only initialize Point<3> objects using the constructor "
425 "that takes three arguments. Point<dim> objects with dim!=3 "
426 "require initialization with the constructor that takes 'dim' "
427 "arguments."));
428
429 // we can only get here if we pass the assertion. use the indirection anyway
430 // so as to avoid compiler warnings about uninitialized elements or writing
431 // beyond the end of the 'values' array
432 constexpr unsigned int y_index = (dim < 2) ? 0 : 1;
433 constexpr unsigned int z_index = (dim < 3) ? 0 : 2;
434 this->values[0] = x;
435 this->values[y_index] = y;
436 this->values[z_index] = z;
437}
438
439
440
441template <int dim, typename Number>
443template <std::size_t dummy_dim,
444 std::enable_if_t<(dim == dummy_dim) && (dummy_dim != 0), int>>
445constexpr Point<dim, Number>::Point(
446 const boost::geometry::model::
447 point<Number, dummy_dim, boost::geometry::cs::cartesian> &boost_pt)
448{
449 Assert(dim <= 3, ExcNotImplemented());
450 this->values[0] = boost::geometry::get<0>(boost_pt);
451 constexpr unsigned int y_index = (dim < 2) ? 0 : 1;
452 constexpr unsigned int z_index = (dim < 3) ? 0 : 2;
453
454 if (dim >= 2)
455 this->values[y_index] = boost::geometry::get<y_index>(boost_pt);
456
457 if (dim >= 3)
458 this->values[z_index] = boost::geometry::get<z_index>(boost_pt);
459}
460
461
462
463template <int dim, typename Number>
465constexpr DEAL_II_HOST_DEVICE
466 Point<dim, Number> Point<dim, Number>::unit_vector(unsigned int i)
467{
469 p[i] = 1.;
470 return p;
471}
472
473
474template <int dim, typename Number>
476constexpr DEAL_II_HOST_DEVICE Number Point<dim, Number>::operator()(
477 const unsigned int index) const
478{
479 AssertIndexRange(static_cast<int>(index), dim);
480 return this->values[index];
481}
482
483
484
485template <int dim, typename Number>
487constexpr DEAL_II_HOST_DEVICE Number &Point<dim, Number>::operator()(
488 const unsigned int index)
489{
490 AssertIndexRange(static_cast<int>(index), dim);
491 return this->values[index];
492}
493
494
495
496template <int dim, typename Number>
498template <typename OtherNumber>
499constexpr DEAL_II_ALWAYS_INLINE Point<dim, Number>
500 &Point<dim, Number>::operator=(const Tensor<1, dim, OtherNumber> &p)
501{
503 return *this;
504}
505
506
507
508template <int dim, typename Number>
510constexpr DEAL_II_HOST_DEVICE Point<dim, Number> Point<dim, Number>::operator+(
511 const Tensor<1, dim, Number> &p) const
512{
513 Point<dim, Number> tmp = *this;
514 tmp += p;
515 return tmp;
516}
517
518
519
520template <int dim, typename Number>
522constexpr DEAL_II_HOST_DEVICE
523 Tensor<1, dim, Number> Point<dim, Number>::operator-(
524 const Point<dim, Number> &p) const
525{
526 return (Tensor<1, dim, Number>(*this) -= p);
527}
528
529
530
531template <int dim, typename Number>
533constexpr DEAL_II_HOST_DEVICE Point<dim, Number> Point<dim, Number>::operator-(
534 const Tensor<1, dim, Number> &p) const
535{
536 Point<dim, Number> tmp = *this;
537 tmp -= p;
538 return tmp;
539}
540
541
542
543template <int dim, typename Number>
545constexpr DEAL_II_HOST_DEVICE Point<dim, Number> Point<dim, Number>::operator-()
546 const
547{
548 Point<dim, Number> result;
549 for (unsigned int i = 0; i < dim; ++i)
550 result.values[i] = -this->values[i];
551 return result;
552}
553
554
555
556template <int dim, typename Number>
558template <typename OtherNumber>
559constexpr DEAL_II_HOST_DEVICE Point<
560 dim,
561 typename ProductType<Number, typename EnableIfScalar<OtherNumber>::type>::
562 type> Point<dim, Number>::operator*(const OtherNumber factor) const
563{
565 for (unsigned int i = 0; i < dim; ++i)
566 tmp[i] = this->operator[](i) * factor;
567 return tmp;
568}
569
570
571
572template <int dim, typename Number>
574template <typename OtherNumber>
575constexpr DEAL_II_HOST_DEVICE Point<
576 dim,
577 typename ProductType<Number, typename EnableIfScalar<OtherNumber>::type>::
578 type> Point<dim, Number>::operator/(const OtherNumber factor) const
579{
580 const Tensor<1, dim, Number> &base_object = *this;
581 return Point<
582 dim,
583 typename ProductType<Number,
584 typename EnableIfScalar<OtherNumber>::type>::type>(
585 ::operator/(base_object, factor));
586}
587
588
589
590template <int dim, typename Number>
592constexpr DEAL_II_HOST_DEVICE Number Point<dim, Number>::operator*(
593 const Tensor<1, dim, Number> &p) const
594{
595 Number res = Number();
596 for (unsigned int i = 0; i < dim; ++i)
597 res += this->operator[](i) * p[i];
598 return res;
599}
600
601
602template <int dim, typename Number>
604constexpr DEAL_II_HOST_DEVICE typename numbers::NumberTraits<Number>::real_type
605 Point<dim, Number>::square() const
606{
607 return this->norm_square();
608}
609
610
611
612template <int dim, typename Number>
614inline DEAL_II_HOST_DEVICE typename numbers::NumberTraits<Number>::real_type
615 Point<dim, Number>::distance(const Point<dim, Number> &p) const
616{
617 return std::sqrt(distance_square(p));
618}
619
620
621
622template <int dim, typename Number>
624constexpr DEAL_II_HOST_DEVICE typename numbers::NumberTraits<Number>::real_type
625 Point<dim, Number>::distance_square(const Point<dim, Number> &p) const
626{
628 for (unsigned int i = 0; i < dim; ++i)
629 {
630 const Number diff = static_cast<Number>(this->values[i]) - p[i];
632 }
633
634 return sum;
635}
636
637
638
639template <int dim, typename Number>
641template <class Archive>
642inline void Point<dim, Number>::serialize(Archive &ar, const unsigned int)
643{
644 // forward to serialization
645 // function in the base class
646 ar &static_cast<Tensor<1, dim, Number> &>(*this);
647}
648
649#endif // DOXYGEN
650
651
652/*--------------------------- Global functions: Point -----------------------*/
653
654
662template <int dim, typename Number, typename OtherNumber>
663constexpr DEAL_II_HOST_DEVICE
664 Point<dim,
665 typename ProductType<Number,
666 typename EnableIfScalar<OtherNumber>::type>::type>
667 operator*(const OtherNumber factor, const Point<dim, Number> &p)
668{
669 return p * factor;
670}
671
672
673
681template <int dim, typename Number>
683inline std::ostream &
684operator<<(std::ostream &out, const Point<dim, Number> &p)
685{
686 for (unsigned int i = 0; i < dim - 1; ++i)
687 out << p[i] << ' ';
688 out << p[dim - 1];
689
690 return out;
691}
692
693
694
701template <int dim, typename Number>
703inline std::istream &
704operator>>(std::istream &in, Point<dim, Number> &p)
705{
706 for (unsigned int i = 0; i < dim; ++i)
707 in >> p[i];
708
709 return in;
710}
711
712
713#ifndef DOXYGEN
714
720template <typename Number>
721inline std::ostream &
722operator<<(std::ostream &out, const Point<1, Number> &p)
723{
724 out << p[0];
725
726 return out;
727}
728
729#endif // DOXYGEN
731
732#endif
std::ostream & operator<<(std::ostream &out, const Point< dim, Number > &p)
Definition point.h:684
Definition point.h:113
constexpr Point< dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const OtherNumber) const
constexpr Point< dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const OtherNumber) const
constexpr Tensor< 1, dim, Number > operator-(const Point< dim, Number > &) const
constexpr Point(const boost::geometry::model::point< Number, dummy_dim, boost::geometry::cs::cartesian > &boost_pt)
constexpr Point< dim, Number > & operator=(const Tensor< 1, dim, OtherNumber > &p)
constexpr numbers::NumberTraits< Number >::real_type distance_square(const Point< dim, Number > &p) const
constexpr Number & operator()(const unsigned int index)
constexpr Point(const Number x, const Number y)
constexpr Point< dim, Number > operator+(const Tensor< 1, dim, Number > &) const
constexpr Point< dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const OtherNumber factor, const Point< dim, Number > &p)
Definition point.h:667
constexpr Point(const Number x, const Number y, const Number z)
static constexpr Point< dim, Number > unit_vector(const unsigned int i)
constexpr Point< dim, Number > operator-() const
constexpr Point(const Number x)
void serialize(Archive &ar, const unsigned int version)
constexpr Point(const Tensor< 1, dim, Number > &)
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
constexpr Point()
constexpr numbers::NumberTraits< Number >::real_type square() const
constexpr Number operator()(const unsigned int index) const
constexpr Point< dim, Number > operator-(const Tensor< 1, dim, Number > &) const
constexpr Number operator*(const Tensor< 1, dim, Number > &p) const
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
Definition tensor.h:852
constexpr Tensor & operator=(const Tensor< rank_, dim, OtherNumber > &rhs)
#define DEAL_II_ALWAYS_INLINE
Definition config.h:161
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:242
#define DEAL_II_HOST_DEVICE
Definition config.h:166
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcMessage(std::string arg1)
T sum(const T &t, const MPI_Comm mpi_communicator)
STL namespace.
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
static constexpr const T & value(const T &t)
Definition numbers.h:683
static constexpr real_type abs_square(const number &x)
Definition numbers.h:565