1 /*
2   HMat-OSS (HMatrix library, open source software)
3 
4   Copyright (C) 2014-2015 Airbus Group SAS
5 
6   This program is free software; you can redistribute it and/or
7   modify it under the terms of the GNU General Public License
8   as published by the Free Software Foundation; either version 2
9   of the License, or (at your option) any later version.
10 
11   This program is distributed in the hope that it will be useful,
12   but WITHOUT ANY WARRANTY; without even the implied warranty of
13   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14   GNU General Public License for more details.
15 
16   You should have received a copy of the GNU General Public License
17   along with this program; if not, write to the Free Software
18   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
19 
20   http://github.com/jeromerobert/hmat-oss
21 */
22 
23 /*! \file
24   \ingroup HMatrix
25   \brief Geometric coordinates.
26 */
27 #ifndef _COORDINATES_HPP
28 #define _COORDINATES_HPP
29 
30 #include <cmath>
31 #include <cstddef>
32 #include <algorithm>
33 
34 namespace hmat {
35 
36 /*! \brief Class representing a 3D point.
37 
38   Used only in examples (cholesky.cpp and kriging.cpp)
39  */
40 class Point {
41 public:
42   union {
43     struct {
44       double x, y, z; /// Coordinates
45     };
46     double xyz[3]; /// Idem
47   };
Point(double _x=0.,double _y=0.,double _z=0.)48   Point(double _x = 0., double _y = 0., double _z = 0.) : x(_x), y(_y), z(_z) {}
49   /*! \brief Return d(this, other)
50    */
distanceTo(const Point & other) const51   inline double distanceTo(const Point &other) const {
52     double result = 0.;
53     double difference = 0.;
54     difference = x - other.x;
55     result += difference * difference;
56     difference = y - other.y;
57     result += difference * difference;
58     difference = z - other.z;
59     result += difference * difference;
60     return sqrt(result);
61   }
62 };
63 
64 
65 class DofCoordinates {
66 public:
67   /*! \brief Create dof coordinates.
68 
69       \param coord  geometrical coordinates
70       \param dim    spatial dimension
71       \param size   number of points
72       \param ownsMemory if true, coordinates are copied
73    */
74   //TODO: why do we always copy (aka ownsMemory is always true)
75   DofCoordinates(double* coord, unsigned dim, unsigned size, bool ownsMemory = false,
76                  unsigned number_of_dof=0, unsigned * span_offsets = NULL, unsigned * spans = NULL);
77 
78   /*! \brief Copy constructor.
79 
80       \param other  instance being copied
81    */
82   DofCoordinates(const DofCoordinates& other);
83 
84   /*! \brief Destructor.
85    */
86   ~DofCoordinates();
87 
88   /*! \brief Get number of points.
89    * \deprecated Legacy API which does not support spans
90    */
91   int size() const;
92 
93   /*! \brief Get spatial dimension.
94    */
dimension() const95   inline int dimension() const { return dimension_; }
96 
97   /*! \brief Accessor to array element.
98    * \deprecated Legacy API which does not support spans
99    */
100   double& get(int i, int j);
101 
102   /*! \brief Accessor to array element.
103    * \deprecated Legacy API which does not support spans
104    */
105   const double& get(int i, int j) const;
106 
numberOfPoints() const107   unsigned numberOfPoints() const { return size_; }
numberOfDof() const108   unsigned numberOfDof() const {
109       return spanOffsets_ == NULL ? size_ : numberOfDof_;
110   }
111 
112   /** Enlarge the given AABB so it include the given DOF */
spanAABB(unsigned dof,double * bb) const113   void spanAABB(unsigned dof, double * bb) const {
114       int n = dof * dimension_;
115       if(spanOffsets_ == NULL) {
116           for (unsigned dim = 0; dim < dimension_; ++dim) {
117               double v = v_[n + dim];
118               bb[dim] = std::min(bb[dim], v);
119               bb[dim + dimension_] = std::max(bb[dim + dimension_], v);
120           }
121       } else {
122           double * aabb = spanAABBs_ + 2 * n;
123           for (unsigned dim = 0; dim < dimension_; ++dim) {
124               bb[dim] = std::min(bb[dim], aabb[dim]);
125               bb[dim + dimension_] = std::max(bb[dim + dimension_],
126                                      aabb[dim + dimension_]);
127           }
128       }
129   }
130 
spanDiameter(unsigned dof,int dim) const131   double spanDiameter(unsigned dof, int dim) const {
132       double d = 0;
133       if(spanOffsets_ != NULL) {
134           double * aabb = spanAABBs_ + 2 * dof * dimension_;
135           d = std::max(aabb[dim + dimension_] - aabb[dim], d);
136       }
137       return d;
138   }
139 
spanSize(unsigned dof) const140   unsigned spanSize(unsigned dof) const {
141       if(spanOffsets_ == NULL)
142           return 1;
143       else
144           return dof == 0 ? spanOffsets_[0] :
145               spanOffsets_[dof] - spanOffsets_[dof - 1];
146   }
147 
spanPoint(unsigned dof,unsigned pointId,unsigned dim) const148   double spanPoint(unsigned dof, unsigned pointId, unsigned dim) const {
149       if(spanOffsets_ == NULL) {
150           return v_[dof * dimension_ + dim];
151       } else {
152           unsigned offset = dof == 0 ? 0 : spanOffsets_[dof - 1];
153           return v_[spans_[offset + pointId] * dimension_ + dim];
154       }
155   }
156 
spanCenter(unsigned dof,unsigned dim) const157   double spanCenter(unsigned dof, unsigned dim) const {
158       if(spanOffsets_ == NULL)
159           return v_[dof * dimension_ + dim];
160       else {
161           double * aabb = spanAABBs_ + dof * dimension_ * 2;
162           return (aabb[dim] + aabb[dim + dimension_]) / 2;
163       }
164   }
165 
166 private:
167   /// Array of size size_ * dimension_
168   double* v_;
169 
170   /// Spatial dimension
171   unsigned dimension_;
172 
173   /// Number of points
174   unsigned size_;
175 
176   /// Flag to tell whether array has been copied and must be freed by destructor
177   const bool ownsMemory_;
178 
179   unsigned numberOfDof_;
180 
181   // Array of size numberOfDof_, giving the position in spans_[] of each dof definition
182   unsigned * spanOffsets_;
183 
184   // Array of size spanOffsets_[numberOfDof_-1], giving the position in v_[] of the span of each dof
185   unsigned * spans_;
186   /**
187    * @brief min_x, miny, ..., maxx, maxy, ... for each DOF.
188    * This is NULL if span is not supportedx else the size is
189    * 2*dimension()*numberOfDof_.
190    */
191   double * spanAABBs_;
192   void init(double* coord, unsigned * span_offsets, unsigned * spans);
193 };
194 
195 } // end namespace hmat
196 
197 #endif
198