diff --git a/Common/OpenFile.cpp b/Common/OpenFile.cpp
index 3f2adb4ccb17d7054117caba82df48437d7328ff..59a50d34a152b6e43028921561cb1c42d3d21ceb 100644
--- a/Common/OpenFile.cpp
+++ b/Common/OpenFile.cpp
@@ -202,7 +202,7 @@ void ParseString(std::string str)
std::string fileName = CTX::instance()->homeDir + CTX::instance()->tmpFileName;
FILE *fp = fopen(fileName.c_str(), "w");
if(fp){
- fprintf(fp, (str + "\n").c_str());
+ fprintf(fp, "%s", (str + "\n").c_str());
fclose(fp);
ParseFile(fileName, true);
GModel::current()->importGEOInternals();
diff --git a/Geo/STensor3.cpp b/Geo/STensor3.cpp
index 216c95517c26149293dee16f8a2518eac0fab005..bafc15d0aa2991ee6096ec7dc547fdb8408733ec 100644
--- a/Geo/STensor3.cpp
+++ b/Geo/STensor3.cpp
@@ -11,12 +11,11 @@ void SMetric3::print (const char *s) const
SMetric3 intersection (const SMetric3 &m1, const SMetric3 &m2)
{
-
SMetric3 im1 = m1.invert();
gmshMatrix<double> V(3,3);
gmshVector<double> S(3);
im1 *= m2;
- im1.eig(V,S);
+ im1.eig(V,S,true);
SVector3 v0(V(0,0),V(1,0),V(2,0));
SVector3 v1(V(0,1),V(1,1),V(2,1));
SVector3 v2(V(0,2),V(1,2),V(2,2));
@@ -57,3 +56,25 @@ SMetric3 interpolation (const SMetric3 &m1,
return im1.invert();
}
+SMetric3 interpolation (const SMetric3 &m1,
+ const SMetric3 &m2,
+ const SMetric3 &m3,
+ const SMetric3 &m4,
+ const double u,
+ const double v,
+ const double w)
+{
+ SMetric3 im1 = m1.invert();
+ SMetric3 im2 = m2.invert();
+ SMetric3 im3 = m3.invert();
+ SMetric3 im4 = m4.invert();
+ im1 *= (1.-u-v-w);
+ im2 *= u;
+ im3 *= v;
+ im4 *= w;
+ im1 += im2;
+ im1 += im3;
+ im1 += im4;
+ return im1.invert();
+}
+
diff --git a/Geo/STensor3.h b/Geo/STensor3.h
index 40b29c703026b4ff39d7effbf89f6b6f09f5eb8b..b7b406b0be5af42e7a4c2cede10760dd7426d9ad 100644
--- a/Geo/STensor3.h
+++ b/Geo/STensor3.h
@@ -32,29 +32,69 @@ class SMetric3 {
for (int j=i;j<3;j++)
_val[getIndex(i,j)] = mat(i,j);
}
+ SMetric3(const SMetric3& m){
+ for (int i=0; i<6; i++) _val[i]=m._val[i];
+ }
// default constructor, identity
SMetric3(const double v = 1.0){
- _val[0] = _val[2] = _val[5] = 1./(v*v);
+ _val[0] = _val[2] = _val[5] = v;
_val[1] = _val[3] = _val[4] = 0.0;
}
- SMetric3(const double l1,
+/* SMetric3(const double l1, // (h_1)^-2 */
+/* const double l2, */
+/* const double l3, */
+/* const SVector3 &t1, */
+/* const SVector3 &t2, */
+/* const SVector3 &t3){ */
+/* SVector3 t1b = t1; */
+/* SVector3 t2b = t2; */
+/* SVector3 t3b = t3; */
+/* t1b[0] *= l1; t2b[0] *= l1; t3b[0] *= l1; */
+/* t1b[1] *= l2; t2b[1] *= l2; t3b[1] *= l2; */
+/* t1b[2] *= l3; t2b[2] *= l3; t3b[2] *= l3; */
+/* _val[0] = dot (t1b,t1); */
+/* _val[1] = dot (t2b,t1); */
+/* _val[2] = dot (t2b,t2); */
+/* _val[3] = dot (t3b,t1); */
+/* _val[4] = dot (t3b,t2); */
+/* _val[5] = dot (t3b,t3); */
+/* } */
+ SMetric3(const double l1, // l1 = h1^-2
const double l2,
const double l3,
const SVector3 &t1,
const SVector3 &t2,
- const SVector3 &t3){
- SVector3 t1b = t1;
- SVector3 t2b = t2;
- SVector3 t3b = t3;
- t1b[0] *= l1; t2b[0] *= l1; t3b[0] *= l1;
- t1b[1] *= l2; t2b[1] *= l2; t3b[1] *= l2;
- t1b[2] *= l3; t2b[2] *= l3; t3b[2] *= l3;
- _val[0] = dot (t1b,t1);
- _val[1] = dot (t2b,t1);
- _val[2] = dot (t2b,t2);
- _val[3] = dot (t3b,t1);
- _val[4] = dot (t3b,t2);
- _val[5] = dot (t3b,t3);
+ const SVector3 &t3)
+ {
+ // M = e^1 * diag * e^1^t
+ // where the elements of diag are l_i = h_i^-2
+ // and the rows of e are the UNIT and ORTHOGONAL directions
+
+ gmshMatrix<double> e(3,3);
+ e(0,0) = t1(0); e(0,1) = t1(1); e(0,2) = t1(2);
+ e(1,0) = t2(0); e(1,1) = t2(1); e(1,2) = t2(2);
+ e(2,0) = t3(0); e(2,1) = t3(1); e(2,2) = t3(2);
+ e.invertInPlace();
+
+ gmshMatrix<double> tmp(3,3);
+ tmp(0,0) = l1 * e(0,0);
+ tmp(0,1) = l2 * e(0,1);
+ tmp(0,2) = l3 * e(0,2);
+ tmp(1,0) = l1 * e(1,0);
+ tmp(1,1) = l2 * e(1,1);
+ tmp(1,2) = l3 * e(1,2);
+ tmp(2,0) = l1 * e(2,0);
+ tmp(2,1) = l2 * e(2,1);
+ tmp(2,2) = l3 * e(2,2);
+
+ e.transposeInPlace();
+
+ _val[0] = tmp(0,0) * e(0,0) + tmp(0,1) * e(1,0) + tmp(0,2) * e(2,0);
+ _val[1] = tmp(1,0) * e(0,0) + tmp(1,1) * e(1,0) + tmp(1,2) * e(2,0);
+ _val[2] = tmp(1,0) * e(0,1) + tmp(1,1) * e(1,1) + tmp(1,2) * e(2,1);
+ _val[3] = tmp(2,0) * e(0,0) + tmp(2,1) * e(1,0) + tmp(2,2) * e(2,0);
+ _val[4] = tmp(2,0) * e(0,1) + tmp(2,1) * e(1,1) + tmp(2,2) * e(2,1);
+ _val[5] = tmp(2,0) * e(0,2) + tmp(2,1) * e(1,2) + tmp(2,2) * e(2,2);
}
inline double &operator()(int i, int j){
return _val[getIndex(i,j)];
@@ -100,11 +140,12 @@ class SMetric3 {
SMetric3 a; a.setMat(result);
return a;
}
- void eig (gmshMatrix<double> &V, gmshVector<double> &S) const {
+ // s: true if eigenvalues are sorted (from min to max of the REAL part)
+ void eig (gmshMatrix<double> &V, gmshVector<double> &S, bool s=false) const {
gmshMatrix<double> me(3,3),right(3,3);
gmshVector<double> im(3);
getMat(me);
- me.eig(V,S,im,right);
+ me.eig(V,S,im,right,s);
}
void print(const char *) const;
};
@@ -112,9 +153,9 @@ class SMetric3 {
// scalar product with respect to the metric
inline double dot(const SVector3 &a, const SMetric3 &m, const SVector3 &b)
{ return
- b.x() * ( m(0,0) * a.x() + m(0,1) * a.y() + m(0,2) * a.z() ) +
- b.y() * ( m(1,0) * a.x() + m(1,1) * a.y() + m(1,2) * a.z() ) +
- b.z() * ( m(2,0) * a.x() + m(2,1) * a.y() + m(2,2) * a.z() ) ;}
+ b.x() * ( m(0,0) * a.x() + m(1,0) * a.y() + m(2,0) * a.z() ) +
+ b.y() * ( m(0,1) * a.x() + m(1,1) * a.y() + m(2,1) * a.z() ) +
+ b.z() * ( m(0,2) * a.x() + m(1,2) * a.y() + m(2,2) * a.z() ) ;}
// compute the largest inscribed ellipsoid...
SMetric3 intersection (const SMetric3 &m1,
@@ -127,5 +168,12 @@ SMetric3 interpolation (const SMetric3 &m1,
const SMetric3 &m3,
const double u,
const double v);
+SMetric3 interpolation (const SMetric3 &m1,
+ const SMetric3 &m2,
+ const SMetric3 &m3,
+ const SMetric3 &m4,
+ const double u,
+ const double v,
+ const double w);
#endif
diff --git a/Numeric/GmshMatrix.cpp b/Numeric/GmshMatrix.cpp
index f8a21ad5ab50dbc9a40ff640198b7d416bf68d67..125933b37287302d7edc17dd1145aad81a6ab6c1 100644
--- a/Numeric/GmshMatrix.cpp
+++ b/Numeric/GmshMatrix.cpp
@@ -143,7 +143,8 @@ template<>
bool gmshMatrix<double>::eig(gmshMatrix<double> &VL, // left eigenvectors
gmshVector<double> &DR, // Real part of eigenvalues
gmshVector<double> &DI, // Im part of eigenvalues
- gmshMatrix<double> &VR )
+ gmshMatrix<double> &VR,
+ bool sortRealPart) // if true: sorted from min '|DR|' to max '|DR|'
{
int N = size1(), info;
int LWORK = 10*N;
@@ -158,12 +159,35 @@ bool gmshMatrix<double>::eig(gmshMatrix<double> &VL, // left eigenvectors
delete [] work;
- if(info == 0) return true;
if(info > 0)
Msg::Error("QR Algorithm failed to compute all the eigenvalues", info, info);
- else
+ else if(info < 0)
Msg::Error("Wrong %d-th argument in eig", -info);
- return false;
+
+ if (sortRealPart) {
+ double tmp[8];
+ // do permutations
+ for (int i=0; i<(size1()-1); i++) {
+ int minR = i;
+ for (int j=i+1; j<size1(); j++) if ( fabs(DR(j)) < fabs(DR(minR)) ) minR = j;
+ if ( minR != i )
+ {
+ tmp[0] = DR(i); tmp[1] = DI(i);
+ tmp[2] = VL(0,i); tmp[3] = VL(1,i); tmp[4] = VL(2,i);
+ tmp[5] = VR(0,i); tmp[6] = VR(1,i); tmp[7] = VR(2,i);
+
+ DR(i) = DR(minR); DI(i) = DI(minR);
+ VL(0,i) = VL(0,minR); VL(1,i) = VL(1,minR); VL(2,i) = VL(2,minR);
+ VR(0,i) = VR(0,minR); VR(1,i) = VR(1,minR); VR(2,i) = VR(2,minR);
+
+ DR(minR) = tmp[0]; DI(minR) = tmp[1];
+ VL(0,minR) = tmp[2]; VL(1,minR) = tmp[3]; VL(2,minR) = tmp[4];
+ VR(0,minR) = tmp[5]; VR(1,minR) = tmp[6]; VR(2,minR) = tmp[7];
+ }
+ }
+ }
+
+ return true;
}
diff --git a/Numeric/GmshMatrix.h b/Numeric/GmshMatrix.h
index 9e6da0f6ffd66e55f0ac2e3f8e721900706fb042..a2d9ca46b2c8f2c809ea6b647fbea41d85b7e293 100644
--- a/Numeric/GmshMatrix.h
+++ b/Numeric/GmshMatrix.h
@@ -60,6 +60,14 @@ class gmshVector
for(int i = 0; i < _r; ++i) s += _data[i]*other._data[i];
return s;
}
+ void print(const char * name="") const {
+ printf("Printing vector %s:\n",name);
+ printf(" ");
+ for(int I = 0; I < size(); I++){
+ printf("%12.5E ",(*this)(I));
+ }
+ printf("\n");
+ }
};
@@ -175,6 +183,19 @@ class gmshMatrix
T(j, i) = (*this)(i, j);
return T;
}
+ inline void transposeInPlace()
+ {
+ if ( size1() != size2() ) {
+ Msg::Error("Not a square matrix (size1: %d, size2: %d)",size1(),size2());
+ }
+ scalar t;
+ for(int i = 0; i < size1(); i++)
+ for(int j = 0; j < i; j++) {
+ t = _data[i + _r * j];
+ _data[i + _r * j] = _data[j + _r * i];
+ _data[j + _r * i] = t;
+ }
+ }
bool lu_solve(const gmshVector<scalar> &rhs, gmshVector<scalar> &result)
#if !defined(HAVE_LAPACK)
{
@@ -194,7 +215,8 @@ class gmshMatrix
bool eig(gmshMatrix<scalar> &VL, // left eigenvectors
gmshVector<double> &DR, // Real part of eigenvalues
gmshVector<double> &DI, // Im part of eigen
- gmshMatrix<scalar> &VR)
+ gmshMatrix<scalar> &VR,
+ bool sortRealPart=false) // if true: sorted from min 'DR' to max 'DR'
#if !defined(HAVE_LAPACK)
{
Msg::Error("Eigenvalue computations requires LAPACK");
@@ -239,6 +261,18 @@ class gmshMatrix
}
#endif
;
+ void print(const char * name="") const {
+ printf("Printing matrix %s:\n",name);
+ int ni = size1();
+ int nj = size2();
+ for(int I = 0; I < ni; I++){
+ printf(" ");
+ for(int J = 0; J < nj; J++){
+ printf("%12.5E ",(*this)(I, J));
+ }
+ printf("\n");
+ }
+ }
};
#endif