← Back to branch summary

~airpollution/fluidity/fluidity_airpollution

« back to all changes in this revision

Viewing changes to libadaptivity/metric_field/include/MetricTensor.h

  • Committer: ziyouzhj
  • Date: 2013-12-09 16:51:29 UTC
  • Revision ID: ziyouzhj@gmail.com-20131209165129-ucoetc3u0atyy05c
airpolution

Show diffs side-by-side

added added

removed removed

Lines of Context:
libadaptivity/metric_field/include/MetricTensor.h
 
1
/* Copyright (C) 2006 Imperial College London and others.
 
2
 
 
3
 Please see the AUTHORS file in the main source directory for a full list
 
4
 of copyright holders.
 
5
 
 
6
 Dr Gerard J Gorman
 
7
 Applied Modelling and Computation Group
 
8
 Department of Earth Science and Engineering
 
9
 Imperial College London
 
10
 
 
11
 g.gorman@imperial.ac.uk
 
12
 
 
13
 This library is free software; you can redistribute it and/or
 
14
 modify it under the terms of the GNU Lesser General Public
 
15
 License as published by the Free Software Foundation; either
 
16
 version 2.1 of the License.
 
17
 
 
18
 This library is distributed in the hope that it will be useful,
 
19
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 
20
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 
21
 Lesser General Public License for more details.
 
22
 
 
23
 You should have received a copy of the GNU Lesser General Public
 
24
 License along with this library; if not, write to the Free Software
 
25
 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
 
26
 USA
 
27
*/
 
28
 
 
29
#ifndef METRICTENSOR_H
 
30
#define METRICTENSOR_H
 
31
 
 
32
#include "confdefs.h"
 
33
 
 
34
#include "vtk.h"
 
35
 
 
36
#include <algorithm>
 
37
#include <cassert>
 
38
#include <cmath>
 
39
#include <iostream>
 
40
#include <string>
 
41
#include <vector>
 
42
 
 
43
#include <float.h>
 
44
 
 
45
#include "cxxblas.h"
 
46
 
 
47
#ifdef FLT_SMALL
 
48
#error FLT_SMALL is already defined
 
49
#endif
 
50
#define FLT_SMALL (1.0E-12)
 
51
 
 
52
/** Performs basic operations on metric tensors.
 
53
 */
 
54
class MetricTensor{
 
55
 public:
 
56
  /// Constructor.
 
57
  MetricTensor(int dimension);
 
58
 
 
59
  /// Constructor.
 
60
  MetricTensor(double m11,
 
61
               double m12, double m22);
 
62
 
 
63
  /// Constructor.
 
64
  MetricTensor(double m11,
 
65
               double m12, double m22,
 
66
               double m13, double m23, double m33);
 
67
 
 
68
  /// Constructor - T stores the bottom triangle of the tensor.
 
69
  MetricTensor(int dim, const double *T);
 
70
 
 
71
  /// Copy constructor.
 
72
  MetricTensor(const MetricTensor&);
 
73
 
 
74
  /// Overloaded assignment operator.
 
75
  const MetricTensor& operator=(const MetricTensor &); 
 
76
 
 
77
  /// Overloaded write operator.
 
78
  friend std::ostream &operator<< (std::ostream&, const MetricTensor&);
 
79
 
 
80
  /// Aspect ratio between the maximum and minimum edge lengths
 
81
  /// required by the metric.
 
82
  double aspect_ratio() const;
 
83
 
 
84
  /// Average desired edge length
 
85
  double average_length() const;
 
86
  
 
87
  /// Calculates metric tensor that circumscribes this tensor and
 
88
  /// input, M. This limits the minimum edge lengths anisotropically.
 
89
  int circumscribe(const MetricTensor& M);
 
90
 
 
91
  /// Determinant of tensor
 
92
  double det() const;
 
93
 
 
94
  /// Calculate eivenvalue decomposition of metric tensor.
 
95
  /// @param D Eigenvalues
 
96
  /// @param V Eigenvectors.
 
97
  int eigen_decomp(double *D, double *V) const;
 
98
 
 
99
  /// Construct metric tensor from eigenvalue decomposition.
 
100
  /// @param D Eigenvalues
 
101
  /// @param V Eigenvectors.
 
102
  int eigen_undecomp(const double *D, const double *V);
 
103
 
 
104
  /// Calculate the inner product between this tensor and M2
 
105
  /// @returns Inner product.
 
106
  MetricTensor dot(const MetricTensor &M2);
 
107
 
 
108
  /// Pass back a copy of the bottom triangle of metric tensor into M (must be allocated by caller).
 
109
  void get_metric(double *M) const;
 
110
 
 
111
  /// Pass back a copy of the full metric tensor into M (must be allocated by caller).
 
112
  void get_metric2(double *M) const;
 
113
 
 
114
  /// Calculates metric tensor that inscribes this tensor and
 
115
  /// input, M. This limits the maximum edge lengths anisotropically.
 
116
  int inscribe(const MetricTensor& M);
 
117
 
 
118
  /// Limit the maximum edge length.
 
119
  int limit_max_size(const double *max_d);
 
120
 
 
121
  /// Limit the minimum edge length.
 
122
  int limit_min_size(const double *min_d);
 
123
 
 
124
  /// Calculate inverse of metric tensor.
 
125
  MetricTensor inv() const;
 
126
 
 
127
  /// True if metric is null.
 
128
  bool is_null() const;
 
129
 
 
130
  /// Limit the aspect ratio of the desired edge lengths.
 
131
  int limit_aspect(double);
 
132
 
 
133
  /// Scale the metric tensor.
 
134
  /// @param s Scale factor
 
135
  void scale(double s);
 
136
 
 
137
  /// Superimposes two metrics.
 
138
  int superimpose(const MetricTensor& M);
 
139
 
 
140
  /// Disable verbose mode.
 
141
  void verbose_off();
 
142
  
 
143
  /// Enable verbose mode.
 
144
  void verbose_on();
 
145
  
 
146
  /// Dump a copy of the metric as vtkPolyData. Must be configured with VTK support.
 
147
  int write_vtk(const char *filename) const;
 
148
 
 
149
 private:
 
150
  size_t dim, t_size;
 
151
  std::vector<double> metric;
 
152
  static bool verbose;
 
153
 
 
154
  /// BLAS [DS]PEV abstraction
 
155
  int blas_spev(char jobz, char uplo, int N, const double ap[],  
 
156
                double eigenvalues[],  double eigenvectors[],  int ldz,  double work[]) const;
 
157
 
 
158
  /// Cofactor of index i, j
 
159
  double cofactor(int i, int j) const;
 
160
 
 
161
  /// Calculate eivenvalue decomposition of tensor T.
 
162
  /// @param T tensor
 
163
  /// @param D Eigenvalues
 
164
  /// @param V Eigenvectors.
 
165
  int eigen_decomp(const double *T, double *D, double *V) const;
 
166
 
 
167
  /// Construct metric tensor from eigenvalue decomposition.
 
168
  /// @param T Tensor
 
169
  /// @param D Eigenvalues
 
170
  /// @param V Eigenvectors.
 
171
  int eigen_undecomp(const double *D, const double *V, double *T) const;
 
172
 
 
173
 
 
174
  /// In-place inversion of tensor.
 
175
  int inv(double&, double&,
 
176
          double&, double&) const;
 
177
 
 
178
  /// In-place inversion of tensor.
 
179
  int inv(double&, double&, double&,
 
180
          double&, double&, double&,
 
181
          double&, double&, double&) const;
 
182
  
 
183
  /// Convert i, j index into lower-triangle index
 
184
  int lookup(size_t i, size_t j) const;
 
185
 
 
186
  /// Map self to the same space where M is the unit sphere
 
187
  int map(const MetricTensor& M);
 
188
 
 
189
  int MatrixDotMatrix(double *A, bool aT, double *B, bool bT, double *C) const;
 
190
  void print_eigen(double *, double*) const;
 
191
 
 
192
  /// Unmap M back from unit metric space.
 
193
  int unmap(const MetricTensor& M);
 
194
};
 
195
 
 
196
#endif