gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
STensor3.cpp
Go to the documentation of this file.
1 // Gmsh - Copyright (C) 1997-2022 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 
6 #include <algorithm>
7 #include "STensor3.h"
8 #include "fullMatrix.h"
9 
10 SMetric3::SMetric3(const double l1, // l1 = h1^-2
11  const double l2, const double l3, const SVector3 &t1,
12  const SVector3 &t2, const SVector3 &t3)
13 {
14  // M = e^t * diag * e
15  // where the elements of diag are l_i = h_i^-2
16  // and the rows of e are the UNIT and ORTHOGONAL directions
17 
18  fullMatrix<double> e(3, 3);
19  e(0, 0) = t1(0);
20  e(0, 1) = t1(1);
21  e(0, 2) = t1(2);
22  e(1, 0) = t2(0);
23  e(1, 1) = t2(1);
24  e(1, 2) = t2(2);
25  e(2, 0) = t3(0);
26  e(2, 1) = t3(1);
27  e(2, 2) = t3(2);
28  e.transposeInPlace();
29 
30  fullMatrix<double> tmp(3, 3);
31  tmp(0, 0) = l1 * e(0, 0);
32  tmp(0, 1) = l2 * e(0, 1);
33  tmp(0, 2) = l3 * e(0, 2);
34  tmp(1, 0) = l1 * e(1, 0);
35  tmp(1, 1) = l2 * e(1, 1);
36  tmp(1, 2) = l3 * e(1, 2);
37  tmp(2, 0) = l1 * e(2, 0);
38  tmp(2, 1) = l2 * e(2, 1);
39  tmp(2, 2) = l3 * e(2, 2);
40 
41  e.transposeInPlace();
42 
43  _val[0] = tmp(0, 0) * e(0, 0) + tmp(0, 1) * e(1, 0) + tmp(0, 2) * e(2, 0);
44  _val[1] = tmp(1, 0) * e(0, 0) + tmp(1, 1) * e(1, 0) + tmp(1, 2) * e(2, 0);
45  _val[2] = tmp(1, 0) * e(0, 1) + tmp(1, 1) * e(1, 1) + tmp(1, 2) * e(2, 1);
46  _val[3] = tmp(2, 0) * e(0, 0) + tmp(2, 1) * e(1, 0) + tmp(2, 2) * e(2, 0);
47  _val[4] = tmp(2, 0) * e(0, 1) + tmp(2, 1) * e(1, 1) + tmp(2, 2) * e(2, 1);
48  _val[5] = tmp(2, 0) * e(0, 2) + tmp(2, 1) * e(1, 2) + tmp(2, 2) * e(2, 2);
49 }
50 
51 void SMetric3::getMat(fullMatrix<double> &mat) const
52 {
53  for(int i = 0; i < 3; i++) {
54  for(int j = 0; j < 3; j++) { mat(i, j) = _val[getIndex(i, j)]; }
55  }
56 }
57 
58 void SMetric3::setMat(const fullMatrix<double> &mat)
59 {
60  for(int i = 0; i < 3; i++)
61  for(int j = 0; j < 3; j++) _val[getIndex(i, j)] = mat(i, j);
62 }
63 
64 SMetric3 SMetric3::invert() const
65 {
66  fullMatrix<double> m(3, 3);
67  getMat(m);
68  m.invertInPlace();
69  SMetric3 ithis;
70  ithis.setMat(m);
71  return ithis;
72 }
73 
74 double SMetric3::determinant() const
75 {
76  fullMatrix<double> m(3, 3);
77  getMat(m);
78  double det = m.determinant();
79  return det;
80 }
81 
82 SMetric3 &SMetric3::operator*=(const SMetric3 &other)
83 {
84  fullMatrix<double> m1(3, 3), m2(3, 3), m3(3, 3);
85  getMat(m1);
86  other.getMat(m2);
87  m1.mult(m2, m3);
88  setMat(m3);
89  return *this;
90 }
91 
92 SMetric3 SMetric3::transform(fullMatrix<double> &V)
93 {
94  fullMatrix<double> m(3, 3);
95  getMat(m);
96  fullMatrix<double> result(3, 3), temp(3, 3);
97  V.transpose().mult(m, temp);
98  temp.mult(V, result);
99  SMetric3 a;
100  a.setMat(result);
101  return a;
102 }
103 
104 void SMetric3::eig(fullMatrix<double> &V, fullVector<double> &S, bool s) const
105 {
106  fullMatrix<double> me(3, 3), left(3, 3);
107  fullVector<double> im(3);
108  getMat(me);
109  me.eig(S, im, left, V, s);
110 }
111 
112 void SMetric3::print(const char *s) const
113 {
114  printf(" metric %s : %12.5E %12.5E %12.5E %12.5E %12.5E %12.5E \n", s,
115  (*this)(0, 0), (*this)(1, 1), (*this)(2, 2), (*this)(0, 1),
116  (*this)(0, 2), (*this)(1, 2));
117 }
118 
119 SMetric3 intersection(const SMetric3 &m1, const SMetric3 &m2)
120 {
121  SMetric3 im1 = m1.invert();
122  fullMatrix<double> V(3, 3);
123  fullVector<double> S(3);
124  im1 *= m2;
125  im1.eig(V, S, true);
126  SVector3 v0(V(0, 0), V(1, 0), V(2, 0));
127  SVector3 v1(V(0, 1), V(1, 1), V(2, 1));
128  SVector3 v2(V(0, 2), V(1, 2), V(2, 2));
129  double l0 = std::max(dot(v0, m1, v0), dot(v0, m2, v0));
130  double l1 = std::max(dot(v1, m1, v1), dot(v1, m2, v1));
131  double l2 = std::max(dot(v2, m1, v2), dot(v2, m2, v2));
132 
133  // Correction from the PhD thesis of Frederic Alauzet p.16
134  // If m2 = alpha*m1, then take the largest metric
135  static const double eps = 1.e-2; // Tolerance to detect triple eigenvalue
136  // (i.e. proportional metrics)
137  const double max_eig = std::max(S(0), std::max(S(1), S(2)));
138  const double min_eig = std::min(S(0), std::min(S(1), S(2)));
139  const double range_eig = fabs((max_eig - min_eig) / max_eig);
140  if(range_eig < eps) return (max_eig >= 1.) ? m2 : m1;
141 
142  SMetric3 iv(l0, l1, l2, v0, v1, v2);
143  return iv;
144 }
145 
146 SMetric3 intersection_alauzet(const SMetric3 &m1, const SMetric3 &m2)
147 {
148  SMetric3 im1 = m1.invert();
149  fullMatrix<double> V(3, 3);
150  fullVector<double> S(3);
151  im1 *= m2;
152  im1.eig(V, S, true);
153  SVector3 v0(V(0, 0), V(1, 0), V(2, 0));
154  SVector3 v1(V(0, 1), V(1, 1), V(2, 1));
155  SVector3 v2(V(0, 2), V(1, 2), V(2, 2));
156  double l0 = std::max(dot(v0, m1, v0), dot(v0, m2, v0));
157  double l1 = std::max(dot(v1, m1, v1), dot(v1, m2, v1));
158  double l2 = std::max(dot(v2, m1, v2), dot(v2, m2, v2));
159  SMetric3 iv(l0, l1, l2, v0, v1, v2);
160  return iv;
161 }
162 
163 // preserve orientation of m1 !!!
164 SMetric3 intersection_conserveM1(const SMetric3 &m1, const SMetric3 &m2)
165 {
166  // we should do
167  // return intersection (m1,m2);
168  fullMatrix<double> V(3, 3);
169  fullVector<double> S(3);
170  m1.eig(V, S, false);
171  SVector3 v0(V(0, 0), V(1, 0), V(2, 0));
172  SVector3 v1(V(0, 1), V(1, 1), V(2, 1));
173  SVector3 v2(V(0, 2), V(1, 2), V(2, 2));
174  double l0 = std::max(dot(v0, m1, v0), dot(v0, m2, v0));
175  double l1 = std::max(dot(v1, m1, v1), dot(v1, m2, v1));
176  double l2 = std::max(dot(v2, m1, v2), dot(v2, m2, v2));
177  SMetric3 iv(l0, l1, l2, v0, v1, v2);
178  return iv;
179 }
180 
181 namespace {
182 
183  double anisoRatio2D(double v0x, double v0y, double v0z, double v1x,
184  double v1y, double v1z, double v2x, double v2y,
185  double v2z, double s0, double s1, double s2)
186  {
187  static const double eps = 1.e-8;
188 
189  double ratio;
190  if(v0x * v0x + v0y * v0y + (v0z - 1.) * (v0z - 1.) <
191  eps) // If 1st eigenvect corresponds to z dir. ...
192  ratio =
193  std::min(fabs(s1 / s2),
194  fabs(s2 / s1)); // ... consider ratio of 2nd and 3rd eigenval.
195  else if(v1x * v1x + v1y * v1y + (v1z - 1.) * (v1z - 1.) <
196  eps) // If 2nd eigenvect corresponds to z dir. ...
197  ratio =
198  std::min(fabs(s0 / s2),
199  fabs(s2 / s0)); // ... consider ratio of 1st and 3rd eigenval.
200  else if(v2x * v2x + v2y * v2y + (v2z - 1.) * (v2z - 1.) <
201  eps) // If 3rd eigenvect corresponds to z dir. ...
202  ratio =
203  std::min(fabs(s0 / s1),
204  fabs(s1 / s0)); // ... consider ratio of 1st and 2nd eigenval.
205  else {
206  printf("Error in anisoRatio2D: z direction not found!\n");
207  ratio = 0.;
208  }
209 
210  return ratio;
211  }
212 
213 } // namespace
214 
215 // preserve orientation of the most anisotropic metric !!!
216 SMetric3 intersection_conserve_mostaniso(const SMetric3 &m1, const SMetric3 &m2)
217 {
218  fullMatrix<double> V1(3, 3);
219  fullVector<double> S1(3);
220  m1.eig(V1, S1, true);
221  double ratio1 =
222  fabs(S1(0) / S1(2)); // Minimum ratio because we take sorted eigenvalues
223  fullMatrix<double> V2(3, 3);
224  fullVector<double> S2(3);
225  m2.eig(V2, S2, true);
226  double ratio2 =
227  fabs(S2(0) / S2(2)); // Minimum ratio because we take sorted eigenvalues
228 
229  if(ratio1 < ratio2)
230  return intersection_conserveM1(m1, m2);
231  else
232  return intersection_conserveM1(m2, m1);
233 }
234 
235 // preserve orientation of the most anisotropic metric in 2D!!!
236 SMetric3 intersection_conserve_mostaniso_2d(const SMetric3 &m1,
237  const SMetric3 &m2)
238 {
239  fullMatrix<double> V1(3, 3);
240  fullVector<double> S1(3);
241  m1.eig(V1, S1, false);
242  double ratio1 =
243  anisoRatio2D(V1(0, 0), V1(1, 0), V1(2, 0), V1(0, 1), V1(1, 1), V1(2, 1),
244  V1(0, 2), V1(1, 2), V1(2, 2), S1(0), S1(1), S1(2));
245 
246  fullMatrix<double> V2(3, 3);
247  fullVector<double> S2(3);
248  m2.eig(V2, S2, false);
249  double ratio2 =
250  anisoRatio2D(V2(0, 0), V2(1, 0), V2(2, 0), V2(0, 1), V2(1, 1), V2(2, 1),
251  V2(0, 2), V2(1, 2), V2(2, 2), S2(0), S2(1), S2(2));
252 
253  if(ratio1 < ratio2)
254  return intersection_conserveM1(m1, m2);
255  else
256  return intersection_conserveM1(m2, m1);
257 }
258 
259 // (1-t) * m1 + t * m2
260 SMetric3 interpolation(const SMetric3 &m1, const SMetric3 &m2, const double t)
261 {
262  SMetric3 im1 = m1.invert();
263  SMetric3 im2 = m2.invert();
264  im1 *= (1. - t);
265  im2 *= t;
266  im1 += im2;
267  return im1.invert();
268 }
269 
270 SMetric3 interpolation(const SMetric3 &m1, const SMetric3 &m2,
271  const SMetric3 &m3, const double u, const double v)
272 {
273  SMetric3 im1 = m1.invert();
274  SMetric3 im2 = m2.invert();
275  SMetric3 im3 = m3.invert();
276  im1 *= (1. - u - v);
277  im2 *= u;
278  im3 *= v;
279  im1 += im2;
280  im1 += im3;
281  return im1.invert();
282 }
283 
284 SMetric3 interpolation(const SMetric3 &m1, const SMetric3 &m2,
285  const SMetric3 &m3, const SMetric3 &m4, const double u,
286  const double v, const double w)
287 {
288  SMetric3 im1 = m1.invert();
289  SMetric3 im2 = m2.invert();
290  SMetric3 im3 = m3.invert();
291  SMetric3 im4 = m4.invert();
292  im1 *= (1. - u - v - w);
293  im2 *= u;
294  im3 *= v;
295  im4 *= w;
296  im1 += im2;
297  im1 += im3;
298  im1 += im4;
299  return im1.invert();
300 }
301 
302 void STensor3::getMat(fullMatrix<double> &mat) const
303 {
304  for(int i = 0; i < 3; i++) {
305  for(int j = 0; j < 3; j++) { mat(i, j) = _val[getIndex(i, j)]; }
306  }
307 }
308 
309 void STensor3::setMat(const fullMatrix<double> &mat)
310 {
311  for(int i = 0; i < 3; i++)
312  for(int j = 0; j < 3; j++) _val[getIndex(i, j)] = mat(i, j);
313 }
314 
315 void STensor3::eig(fullMatrix<double> &V, fullVector<double> &S, bool s) const
316 {
317  fullMatrix<double> me(3, 3), left(3, 3);
318  fullVector<double> im(3);
319  this->getMat(me);
320  me.eig(S, im, left, V, s);
321 }
322 
323 void STensor3::print(const char *s) const
324 {
325  printf(
326  " tensor %s : \n"
327  " %12.5E %12.5E %12.5E \n %12.5E %12.5E %12.5E \n %12.5E %12.5E %12.5E \n",
328  s, (*this)(0, 0), (*this)(0, 1), (*this)(0, 2), (*this)(1, 0),
329  (*this)(1, 1), (*this)(1, 2), (*this)(2, 0), (*this)(2, 1), (*this)(2, 2));
330 }
STensor3::getIndex
int getIndex(int i, int j) const
Definition: STensor3.h:114
STensor3::print
void print(const char *) const
Definition: STensor3.cpp:323
dot
double dot(const SVector3 &a, const SMetric3 &m, const SVector3 &b)
Definition: STensor3.h:71
SMetric3
Definition: STensor3.h:17
SMetric3::SMetric3
SMetric3(const double v=1.0)
Definition: STensor3.h:25
fullVector< double >
STensor3::_val
double _val[9]
Definition: STensor3.h:100
SMetric3::eig
void eig(fullMatrix< double > &V, fullVector< double > &S, bool s=false) const
Definition: STensor3.cpp:104
SVector3
Definition: SVector3.h:16
fullMatrix::transpose
fullMatrix< scalar > transpose() const
Definition: fullMatrix.h:558
intersection
SMetric3 intersection(const SMetric3 &m1, const SMetric3 &m2)
Definition: STensor3.cpp:119
SMetric3::invert
SMetric3 invert() const
Definition: STensor3.cpp:64
SMetric3::determinant
double determinant() const
Definition: STensor3.cpp:74
fullMatrix< double >
SMetric3::getMat
void getMat(fullMatrix< double > &mat) const
Definition: STensor3.cpp:51
SMetric3::print
void print(const char *) const
Definition: STensor3.cpp:112
SMetric3::getIndex
int getIndex(int i, int j) const
Definition: STensor3.h:37
SMetric3::_val
double _val[6]
Definition: STensor3.h:21
SMetric3::operator*=
SMetric3 & operator*=(const double &other)
Definition: STensor3.h:59
fullMatrix::invertInPlace
bool invertInPlace()
Definition: fullMatrix.h:662
interpolation
SMetric3 interpolation(const SMetric3 &m1, const SMetric3 &m2, const double t)
Definition: STensor3.cpp:260
intersection_conserve_mostaniso
SMetric3 intersection_conserve_mostaniso(const SMetric3 &m1, const SMetric3 &m2)
Definition: STensor3.cpp:216
fullMatrix::determinant
scalar determinant() const
Definition: fullMatrix.h:681
S
#define S
Definition: DefaultOptions.h:21
fullMatrix::transposeInPlace
void transposeInPlace()
Definition: fullMatrix.h:565
intersection_conserve_mostaniso_2d
SMetric3 intersection_conserve_mostaniso_2d(const SMetric3 &m1, const SMetric3 &m2)
Definition: STensor3.cpp:236
STensor3.h
STensor3::eig
void eig(fullMatrix< double > &V, fullVector< double > &S, bool s=false) const
Definition: STensor3.cpp:315
fullMatrix::mult
void mult(const fullVector< scalar > &x, fullVector< scalar > &y) const
Definition: fullMatrix.h:487
SMetric3::transform
SMetric3 transform(fullMatrix< double > &V)
Definition: STensor3.cpp:92
intersection_alauzet
SMetric3 intersection_alauzet(const SMetric3 &m1, const SMetric3 &m2)
Definition: STensor3.cpp:146
STensor3::getMat
void getMat(fullMatrix< double > &mat) const
Definition: STensor3.cpp:302
intersection_conserveM1
SMetric3 intersection_conserveM1(const SMetric3 &m1, const SMetric3 &m2)
Definition: STensor3.cpp:164
STensor3::setMat
void setMat(const fullMatrix< double > &mat)
Definition: STensor3.cpp:309
SMetric3::setMat
void setMat(const fullMatrix< double > &mat)
Definition: STensor3.cpp:58
fullMatrix.h