diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index 2fa6d434ee055ed2314cef97fd72f1e4d8ca9838..5d92a24bef0187c8208c0c2a45d7790213965912 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -1,6 +1,24 @@
#ifndef _GMSH_DEFINES_H_
#define _GMSH_DEFINES_H_
+// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 // USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
// IO file formats
#define FORMAT_MSH 1
#define FORMAT_UNV 2
diff --git a/Common/GmshMatrix.h b/Common/GmshMatrix.h
index eaa638400192ab9bc9275fb9ac65f0539d44bae2..41a2c069b0ab983040ba1a6d99390375855ef80f 100644
--- a/Common/GmshMatrix.h
+++ b/Common/GmshMatrix.h
@@ -1,5 +1,23 @@
-#ifndef _GMSH_BOOSTMATRIX_
-#define _GMSH_BOOSTMATRIX_
+#ifndef _GMSH_MATRIX_H_
+#define _GMSH_MATRIX_H_
+
+// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 // USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
#include <assert.h>
@@ -9,20 +27,20 @@ class Gmsh_Vector
private:
int r;
public:
- inline int size() const {return r;}
+ inline int size() const { return r; }
SCALAR *data;
- ~Gmsh_Vector() {delete [] data;}
+ ~Gmsh_Vector() { delete [] data; }
Gmsh_Vector(int R)
: r(R)
{
- data = new SCALAR [r];
+ data = new SCALAR[r];
scal(0);
}
Gmsh_Vector(const Gmsh_Vector<SCALAR> &other)
: r(other.r)
{
- data = new double [r];
- copy ( other.data ) ;
+ data = new double[r];
+ copy(other.data);
}
inline SCALAR operator () (int i) const
{
@@ -32,19 +50,19 @@ public:
{
return data[i];
}
- inline SCALAR operator *(const Gmsh_Vector<SCALAR> & other)
+ inline SCALAR operator *(const Gmsh_Vector<SCALAR> &other)
{
throw;
}
- inline void scal ( const SCALAR s)
+ inline void scal(const SCALAR s)
{
- for (int i=0;i<r;++i)data[i]*=s;
+ for (int i = 0; i < r; ++i) data[i] *= s;
}
- inline void copy ( const SCALAR **other)
+ inline void copy(const SCALAR **other)
{
- for (int i=0;i<r;++i)data[i]=other.data[i];
+ for (int i = 0; i < r; ++i) data[i] = other.data[i];
}
- inline void lu_solve (const Gmsh_Vector<SCALAR>& rhs, Gmsh_Vector<SCALAR> & result)
+ inline void lu_solve (const Gmsh_Vector<SCALAR> &rhs, Gmsh_Vector<SCALAR> &result)
{
throw;
}
@@ -54,223 +72,211 @@ template <class SCALAR>
class Gmsh_Matrix
{
private:
- int r,c;
+ int r, c;
public:
- inline int size1() const {return r;}
- inline int size2() const {return c;}
+ inline int size1() const { return r; }
+ inline int size2() const { return c; }
SCALAR *data;
- ~Gmsh_Matrix() {delete [] data;}
+ ~Gmsh_Matrix() { delete [] data; }
Gmsh_Matrix(int R,int C)
- : r(R),c(C)
+ : r(R), c(C)
{
- data = new SCALAR [r*c];
+ data = new SCALAR[r * c];
scal(0);
}
Gmsh_Matrix(const Gmsh_Matrix<SCALAR> &other)
- : r(other.r),c(other.c)
+ : r(other.r), c(other.c)
{
- data = new double [r*c];
- copy ( other.data ) ;
+ data = new double[r * c];
+ copy(other.data);
}
inline SCALAR operator () (int i, int j) const
{
- return data[i+r*j];
+ return data[i + r * j];
}
inline SCALAR & operator () (int i, int j)
{
- return data[i+r*j];
+ return data[i + r * j];
}
- inline Gmsh_Matrix operator *(const Gmsh_Matrix<SCALAR> & other)
+ inline Gmsh_Matrix operator *(const Gmsh_Matrix<SCALAR> &other)
{
throw;
}
- inline void scal ( const SCALAR s)
+ inline void scal(const SCALAR s)
{
- for (int i=0;i<r*c;++i)data[i]*=s;
+ for (int i = 0; i < r * c; ++i) data[i] *= s;
}
- inline void copy ( const SCALAR **other)
+ inline void copy(const SCALAR **other)
{
- for (int i=0;i<r*c;++i)data[i]=other.data[i];
+ for (int i = 0; i < r * c; ++i) data[i] = other.data[i];
}
- inline void mult(const Gmsh_Matrix<SCALAR> & x, const Gmsh_Matrix<SCALAR> & b)
+ inline void mult(const Gmsh_Matrix<SCALAR> &x, const Gmsh_Matrix<SCALAR> &b)
{
throw;
}
- inline void mult (const Gmsh_Vector<SCALAR> & x, Gmsh_Vector<SCALAR> & b )
+ inline void mult (const Gmsh_Vector<SCALAR> &x, Gmsh_Vector<SCALAR> &b)
{
throw;
}
- inline void least_squares (const Gmsh_Vector<SCALAR> & rhs, Gmsh_Vector<SCALAR> & result)
+ inline void least_squares(const Gmsh_Vector<SCALAR> &rhs, Gmsh_Vector<SCALAR> &result)
{
throw;
}
- inline void lu_solve (const Gmsh_Vector<SCALAR> & rhs, Gmsh_Vector<SCALAR> & result)
+ inline void lu_solve (const Gmsh_Vector<SCALAR> &rhs, Gmsh_Vector<SCALAR> &result)
{
throw;
}
- Gmsh_Matrix cofactor(int i,int j) const
+ Gmsh_Matrix cofactor(int i, int j) const
{
throw;
}
- inline void invert ()
+ inline void invert()
{
throw;
}
- double determinant() const {
+ double determinant() const
+ {
throw;
}
};
#ifdef HAVE_GSL
+
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
+
class GSL_Vector
{
private:
int r;
public:
- inline int size() const {return r;}
+ inline int size() const { return r; }
gsl_vector *data;
- ~GSL_Vector() {gsl_vector_free (data);}
- GSL_Vector(int R)
- : r(R)
+ ~GSL_Vector() { gsl_vector_free(data); }
+ GSL_Vector(int R) : r(R)
{
- data = gsl_vector_calloc (r);
+ data = gsl_vector_calloc(r);
}
- GSL_Vector(const GSL_Vector&other)
- : r(other.r)
+ GSL_Vector(const GSL_Vector &other) : r(other.r)
{
- data = gsl_vector_calloc (r);
- gsl_vector_memcpy (data, other.data);
+ data = gsl_vector_calloc(r);
+ gsl_vector_memcpy(data, other.data);
}
inline double operator () (int i) const
{
- return gsl_vector_get (data,i);
+ return gsl_vector_get(data, i);
}
inline double & operator () (int i)
{
- return *gsl_vector_ptr (data,i);
+ return *gsl_vector_ptr(data, i);
}
inline double norm(){
return gsl_blas_dnrm2(data);
}
- inline void scale (const double & y)
- {
- if (y == 0.0) gsl_vector_set_zero( data );
- else gsl_vector_scale ( data, y );
- }
+ inline void scale(const double &y)
+ {
+ if (y == 0.0) gsl_vector_set_zero(data);
+ else gsl_vector_scale(data, y);
+ }
};
class GSL_Matrix
{
private:
gsl_matrix_view view;
- int r,c;
+ int r, c;
public:
- inline size_t size1() const {return data->size1;}
- inline size_t size2() const {return data->size2;}
+ inline size_t size1() const { return data->size1; }
+ inline size_t size2() const { return data->size2; }
gsl_matrix *data;
- GSL_Matrix(gsl_matrix_view _data)
- : view(_data),data(&view.matrix)
- {}
-
- GSL_Matrix(size_t R, size_t C)
- { data = gsl_matrix_calloc (R, C); }
-
- GSL_Matrix() : data(0)
- {}
-
- GSL_Matrix(const GSL_Matrix&other):
- data(0)
- {
- if(data)gsl_matrix_free (data);
- data = gsl_matrix_calloc (other.data->size1, other.data->size2);
- gsl_matrix_memcpy (data, other.data);
- }
- virtual ~GSL_Matrix() {if(data && data->owner == 1)gsl_matrix_free (data);}
+ GSL_Matrix(gsl_matrix_view _data) : view(_data),data(&view.matrix) {}
+ GSL_Matrix(size_t R, size_t C) { data = gsl_matrix_calloc (R, C); }
+ GSL_Matrix() : data(0) {}
+ GSL_Matrix(const GSL_Matrix&other) : data(0)
+ {
+ if(data) gsl_matrix_free(data);
+ data = gsl_matrix_calloc(other.data->size1, other.data->size2);
+ gsl_matrix_memcpy(data, other.data);
+ }
+ virtual ~GSL_Matrix() { if(data && data->owner == 1) gsl_matrix_free(data); }
GSL_Matrix & operator = (const GSL_Matrix&other)
- {
- if (&other != this)
- {
- if(data)gsl_matrix_free (data);
- data = gsl_matrix_calloc (other.data->size1, other.data->size2);
- gsl_matrix_memcpy (data, other.data);
- }
- return *this;
+ {
+ if (&other != this){
+ if(data) gsl_matrix_free(data);
+ data = gsl_matrix_calloc(other.data->size1, other.data->size2);
+ gsl_matrix_memcpy(data, other.data);
}
-
- void memcpy (const GSL_Matrix&other)
+ return *this;
+ }
+ void memcpy(const GSL_Matrix &other)
{
- gsl_matrix_memcpy (data, other.data);
+ gsl_matrix_memcpy(data, other.data);
}
-
inline double operator () (int i, int j) const
{
- return gsl_matrix_get (data,i,j);
+ return gsl_matrix_get(data, i, j);
}
inline double & operator () (int i, int j)
{
- return *gsl_matrix_ptr (data,i,j);
+ return *gsl_matrix_ptr(data, i, j);
}
- inline void mult(const GSL_Matrix & x, const GSL_Matrix & b)
+ inline void mult(const GSL_Matrix &x, const GSL_Matrix &b)
{
- gsl_blas_dgemm (CblasNoTrans,CblasNoTrans, 1.0, data, x.data, 1.0, b.data);
+ gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, data, x.data, 1.0, b.data);
}
- inline void set_all (const double & m )
- {
- gsl_matrix_set_all (data, m);
- }
- inline void least_squares (const GSL_Vector & rhs, GSL_Vector & result)
+ inline void set_all(const double &m)
+ {
+ gsl_matrix_set_all(data, m);
+ }
+ inline void least_squares (const GSL_Vector &rhs, GSL_Vector &result)
{
- assert (r > c);
- assert (rhs.size() == r);
- assert (result.size() == c);
- GSL_Matrix *ls = new GSL_Matrix(c, c);
+ assert(r > c);
+ assert(rhs.size() == r);
+ assert(result.size() == c);
+ GSL_Matrix *ls = new GSL_Matrix(c, c);
GSL_Vector *ls_rhs = new GSL_Vector(c);
- //GSL_Vector *test = new GSL_Vector(c);
- gsl_blas_dgemm (CblasTrans,CblasNoTrans, 1.0, data, data, 1.0, ls->data);
- gsl_blas_dgemv (CblasTrans, 1.0, data, rhs.data, 1.0, ls_rhs->data);
+ gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, data, data, 1.0, ls->data);
+ gsl_blas_dgemv(CblasTrans, 1.0, data, rhs.data, 1.0, ls_rhs->data);
ls->lu_solve (*ls_rhs,result);
delete ls;
delete ls_rhs;
}
- inline void lu_solve (const GSL_Vector & rhs, GSL_Vector & result)
+ inline void lu_solve(const GSL_Vector &rhs, GSL_Vector &result)
{
int s;
- gsl_permutation * p = gsl_permutation_alloc (size1());
- gsl_linalg_LU_decomp ( data, p, &s);
- gsl_linalg_LU_solve ( data , p, rhs.data, result.data ) ;
- gsl_permutation_free (p);
+ gsl_permutation * p = gsl_permutation_alloc(size1());
+ gsl_linalg_LU_decomp(data, p, &s);
+ gsl_linalg_LU_solve(data, p, rhs.data, result.data);
+ gsl_permutation_free(p);
}
inline void invert ()
{
int s;
- gsl_permutation * p = gsl_permutation_alloc (size1());
- gsl_linalg_LU_decomp ( data, p, &s);
- gsl_matrix *data_inv = gsl_matrix_calloc (size1(), size2());
- gsl_linalg_LU_invert ( data , p, data_inv ) ;
- gsl_matrix_memcpy (data, data_inv);
- gsl_matrix_free (data_inv);
+ gsl_permutation *p = gsl_permutation_alloc (size1());
+ gsl_linalg_LU_decomp(data, p, &s);
+ gsl_matrix *data_inv = gsl_matrix_calloc(size1(), size2());
+ gsl_linalg_LU_invert(data, p, data_inv) ;
+ gsl_matrix_memcpy(data, data_inv);
+ gsl_matrix_free(data_inv);
gsl_permutation_free (p);
}
- inline bool invertSecure (double& det)
- {
- int s;
- gsl_permutation * p = gsl_permutation_alloc (size1());
- gsl_linalg_LU_decomp ( data, p, &s);
- det = gsl_linalg_LU_det(data,s);
- gsl_matrix *data_inv = gsl_matrix_calloc (size1(), size2());
- gsl_linalg_LU_invert ( data , p, data_inv ) ;
- gsl_matrix_memcpy (data, data_inv);
- gsl_matrix_free (data_inv);
- gsl_permutation_free (p);
-
- return (det != 0.);
- }
+ inline bool invertSecure(double &det)
+ {
+ int s;
+ gsl_permutation *p = gsl_permutation_alloc (size1());
+ gsl_linalg_LU_decomp(data, p, &s);
+ det = gsl_linalg_LU_det(data, s);
+ gsl_matrix *data_inv = gsl_matrix_calloc(size1(), size2());
+ gsl_linalg_LU_invert(data, p, data_inv);
+ gsl_matrix_memcpy(data, data_inv);
+ gsl_matrix_free(data_inv);
+ gsl_permutation_free(p);
+ return (det != 0.);
+ }
double determinant() const
{
GSL_Matrix copy = *this;
@@ -278,49 +284,54 @@ public:
copy.invertSecure(det);
return det;
}
- GSL_Matrix cofactor(int i,int j) const
- {
- int ni = size1();
- int nj = size2();
- GSL_Matrix cof(ni-1,nj-1);
- if (i>0) {
- if (j>0) GSL_Matrix(cof.touchSubmatrix(0,i,0, j )).memcpy(GSL_Matrix(seeSubmatrix(0,i,0 , j)));
- if (j<nj-1) GSL_Matrix(cof.touchSubmatrix(0,i,j,nj-j-1)).memcpy(GSL_Matrix(seeSubmatrix(0,i,j+1,nj-j-1)));
- }
-
- if (i<ni-1) {
- if (j<nj-1) GSL_Matrix(cof.touchSubmatrix(i,ni-i-1,j,nj-j-1)).memcpy(GSL_Matrix(seeSubmatrix(i+1,ni-i-1,j+1,nj-j-1)));
- if (j>0) GSL_Matrix(cof.touchSubmatrix(i,ni-i-1,0, j)).memcpy(GSL_Matrix(seeSubmatrix(i+1,ni-i-1,0 , j)));
- }
- return cof;
+ GSL_Matrix cofactor(int i, int j) const
+ {
+ int ni = size1();
+ int nj = size2();
+ GSL_Matrix cof(ni - 1, nj - 1);
+ if (i > 0) {
+ if (j > 0)
+ GSL_Matrix(cof.touchSubmatrix(0, i , 0, j)).
+ memcpy(GSL_Matrix(seeSubmatrix(0, i, 0, j)));
+ if (j < nj - 1)
+ GSL_Matrix(cof.touchSubmatrix(0, i, j, nj - j - 1)).
+ memcpy(GSL_Matrix(seeSubmatrix(0, i, j + 1,nj - j - 1)));
}
-
- inline void mult (const GSL_Vector & x, GSL_Vector & b )
+ if (i < ni - 1) {
+ if (j < nj - 1)
+ GSL_Matrix(cof.touchSubmatrix(i, ni - i - 1, j, nj - j - 1)).
+ memcpy(GSL_Matrix(seeSubmatrix(i + 1, ni - i - 1, j + 1, nj - j - 1)));
+ if (j > 0)
+ GSL_Matrix(cof.touchSubmatrix(i, ni - i - 1, 0, j)).
+ memcpy(GSL_Matrix(seeSubmatrix(i + 1, ni - i - 1, 0, j)));
+ }
+ return cof;
+ }
+ inline void mult(const GSL_Vector &x, GSL_Vector &b)
{
- gsl_blas_dgemv (CblasNoTrans, 1.0, data, x.data, 1.0, b.data);
+ gsl_blas_dgemv(CblasNoTrans, 1.0, data, x.data, 1.0, b.data);
}
- inline gsl_matrix_view touchSubmatrix(int i0,int ni,int j0,int nj)
+ inline gsl_matrix_view touchSubmatrix(int i0, int ni, int j0, int nj)
{
- return gsl_matrix_submatrix (data,i0,j0,ni,nj);
+ return gsl_matrix_submatrix(data, i0, j0, ni, nj);
}
- inline const gsl_matrix_view seeSubmatrix(int i0,int ni,int j0,int nj) const
+ inline const gsl_matrix_view seeSubmatrix(int i0, int ni, int j0, int nj) const
{
- return gsl_matrix_submatrix (data,i0,j0,ni,nj);
+ return gsl_matrix_submatrix(data, i0, j0, ni, nj);
+ }
+ inline void scale(const double &m)
+ {
+ if (m == 0.0) gsl_matrix_set_zero(data);
+ else gsl_matrix_scale(data, m);
+ }
+ inline void add(const double &a)
+ {
+ gsl_matrix_add_constant(data, a);
+ }
+ inline void add(const GSL_Matrix &m)
+ {
+ gsl_matrix_add (data, m.data);
}
- inline void scale (const double & m )
- {
- if (m == 0.0) gsl_matrix_set_zero ( data );
- else gsl_matrix_scale (data, m);
- }
- inline void add (const double & a )
- {
- gsl_matrix_add_constant (data, a);
- }
- inline void add (const GSL_Matrix & m )
- {
- gsl_matrix_add (data, m.data);
- }
-
};
typedef GSL_Matrix Double_Matrix;
typedef GSL_Vector Double_Vector;
@@ -328,6 +339,6 @@ typedef GSL_Vector Double_Vector;
typedef Gmsh_Matrix<double> Double_Matrix;
typedef Gmsh_Vector<double> Double_Vector;
#endif
-typedef Gmsh_Matrix<int> Int_Matrix;
-typedef Gmsh_Vector<int> Int_Vector;
+typedef Gmsh_Matrix<int> Int_Matrix;
+typedef Gmsh_Vector<int> Int_Vector;
#endif
diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index 7f035b612ed3c12eac19e273c3fedf8c6a984330..5aad41505182da053f9520d4a7e09033389f2f67 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: GEdge.cpp,v 1.39 2008-02-21 09:45:15 remacle Exp $
+// $Id: GEdge.cpp,v 1.40 2008-02-21 12:11:12 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -25,12 +25,12 @@
#include "GFace.h"
#include "MElement.h"
#include "GmshDefines.h"
-#include "GaussLegendre1D.h"
#if defined(HAVE_GMSH_EMBEDDED)
# include "GmshEmbedded.h"
#else
# include "Message.h"
+# include "GaussLegendre1D.h"
#endif
GEdge::GEdge(GModel *model, int tag, GVertex *_v0, GVertex *_v1)
@@ -172,17 +172,21 @@ double GEdge::curvature(double par) const
return norm(d);
}
-
-double GEdge::length (const double &u0, const double &u1, const int nbQuadPoints){
- double *t=0,*w=0;
- gmshGaussLegendre1D (nbQuadPoints , &t, &w);
+double GEdge::length(const double &u0, const double &u1, const int nbQuadPoints)
+{
+#if defined(HAVE_GMSH_EMBEDDED)
+ return -1.;
+#else
+ double *t = 0, *w = 0;
+ gmshGaussLegendre1D(nbQuadPoints, &t, &w);
double L = 0.0;
- const double rapJ = (u1-u0)*.5;
- for (int i=0;i<nbQuadPoints;i++){
- const double ui = u0 * 0.5 * (1.-t[i]) + u1 * 0.5 * (1.+t[i]);
+ const double rapJ = (u1 - u0) * .5;
+ for (int i = 0; i < nbQuadPoints; i++){
+ const double ui = u0 * 0.5 * (1. - t[i]) + u1 * 0.5 * (1. + t[i]);
SVector3 der = firstDer(ui);
const double d = norm(der);
L += d * w[i] * rapJ;
}
return L;
+#endif
}
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index b52403652c8e8e7e2d5e9418d4b6061d23b76383..487bb1591b805465d783c77b8e3a43d7269f3281 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -102,7 +102,7 @@ class GEdge : public GEntity {
// the length of the model edge
inline double length() const { return _length; }
inline void setLength(const double l) { _length = l; }
- double length (const double &u0, const double &u1, const int nbQuadPoints = 4);
+ double length(const double &u0, const double &u1, const int nbQuadPoints = 4);
// one can impose the mesh size at an edge
virtual double prescribedMeshSizeAtVertex() const { return meshAttributes.meshSize; }
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index dea69e09e5409170a13507df1d32bf1bc843e38a..8aa5c49c7354706f15dfabc458c18382edbb5efd 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1,4 +1,4 @@
-// $Id: GFace.cpp,v 1.52 2008-02-21 09:45:15 remacle Exp $
+// $Id: GFace.cpp,v 1.53 2008-02-21 12:11:12 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -23,13 +23,13 @@
#include "GFace.h"
#include "GEdge.h"
#include "MElement.h"
-#include "GaussLegendre1D.h"
#if defined(HAVE_GMSH_EMBEDDED)
# include "GmshEmbedded.h"
#else
# include "Message.h"
# include "Numeric.h"
+# include "GaussLegendre1D.h"
# include "VertexArray.h"
# include "Context.h"
# if defined(HAVE_GSL)
@@ -485,8 +485,7 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z,
GPoint P = point(U,V);
err2 = sqrt(DSQR(X - P.x()) + DSQR(Y - P.y()) + DSQR(Z - P.z()));
- if (err2 < 1.e-8 * CTX.lc)return;
-
+ if (err2 < 1.e-8 * CTX.lc) return;
while(err > Precision && iter < MaxIter) {
P = point(U,V);
@@ -688,28 +687,34 @@ bool GFace::buildSTLTriangulation()
return true;
#endif
}
+
// by default we assume that straight lines are geodesics
-SPoint2 GFace::geodesic(const SPoint2 &pt1 , const SPoint2 &pt2 , double t){
- return pt1 + (pt2-pt1) * t;
+SPoint2 GFace::geodesic(const SPoint2 &pt1 , const SPoint2 &pt2 , double t)
+{
+ return pt1 + (pt2 - pt1) * t;
}
+
// length of a curve drawn on a surface
// S = (X(u,v), Y(u,v), Z(u,v) );
// u = u(t) , v = v(t)
// C = C ( u(t), v(t) )
// dC/dt = dC/du du/dt + dC/dv dv/dt
-double GFace::length(const SPoint2 &pt1 , const SPoint2 &pt2, int nbQuadPoints){
- double *t=0,*w=0;
+double GFace::length(const SPoint2 &pt1, const SPoint2 &pt2, int nbQuadPoints)
+{
+#if defined(HAVE_GMSH_EMBEDDED)
+ return -1.;
+#else
+ double *t = 0, *w = 0;
double L = 0.0;
- gmshGaussLegendre1D (nbQuadPoints , &t, &w);
- for (int i=0;i<nbQuadPoints;i++){
- const double ti = 0.5 * (1.+t[i]);
- SPoint2 pi = geodesic (pt1, pt2, ti);
- Pair<SVector3,SVector3> der2 = firstDer(pi);
- SVector3 der = der2.left() * (pt2.x()-pt1.x()) + der2.right() * (pt2.y()-pt1.y());
+ gmshGaussLegendre1D(nbQuadPoints, &t, &w);
+ for (int i = 0; i < nbQuadPoints; i++){
+ const double ti = 0.5 * (1. + t[i]);
+ SPoint2 pi = geodesic(pt1, pt2, ti);
+ Pair<SVector3, SVector3> der2 = firstDer(pi);
+ SVector3 der = der2.left() * (pt2.x() - pt1.x()) + der2.right() * (pt2.y() - pt1.y());
const double d = norm(der);
L += d * w[i] ;
}
return L;
-
-
+#endif
}
diff --git a/Geo/GFace.h b/Geo/GFace.h
index fb416da9d1ffc81178e665cad162e4c576e74111..619b56afec097a31fd3213746ce609d5b332ea16 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -100,11 +100,13 @@ class GFace : public GEntity
// Return the point on the face corresponding to the given parameter.
virtual GPoint point(double par1, double par2) const = 0;
virtual GPoint point(const SPoint2 &pt) const { return point(pt.x(), pt.y()); }
+
// computes, in parametric space, the interpolation from pt1 to pt2 alon the geodesic
// of the surface
- virtual SPoint2 geodesic ( const SPoint2 &pt1 , const SPoint2 &pt2 , double t);
+ virtual SPoint2 geodesic(const SPoint2 &pt1, const SPoint2 &pt2, double t);
+
// computes the length of a geodesic between two points in parametric space
- virtual double length ( const SPoint2 &pt1 , const SPoint2 &pt2, int n=4);
+ virtual double length(const SPoint2 &pt1, const SPoint2 &pt2, int n=4);
// If the mapping is a conforming mapping, i.e. a mapping that
// conserves angles, this function returns the eigenvalue of the
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 589d788d3571b7030827baaec85e66229b8b16d2..56f00755ffae010f9997d0ee326f5d710b00f120 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.53 2008-02-21 09:45:15 remacle Exp $
+// $Id: MElement.cpp,v 1.54 2008-02-21 12:11:12 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -23,13 +23,13 @@
#include "MElement.h"
#include "GEntity.h"
#include "GFace.h"
-#include "FunctionSpace.h"
#include "GmshMatrix.h"
#if defined(HAVE_GMSH_EMBEDDED)
# include "GmshEmbedded.h"
#else
# include "Numeric.h"
+# include "FunctionSpace.h"
# include "Message.h"
# include "Context.h"
# include "qualityMeasures.h"
@@ -643,89 +643,91 @@ void MTriangle::circumcenterXY(double *res) const
circumcenterXY(p1, p2, p3, res);
}
-
void MTriangle::jac(int ord, MVertex *vs[], double uu, double vv, double j[2][3])
{
+#if defined(HAVE_GMSH_EMBEDDED)
+ return -1.;
+#else
double grads[256][2];
int nf = getNumFaceVertices();
if (!nf){
switch(ord){
- case 1: gmshFunctionSpaces::find(MSH_TRI_3).df(uu,vv,grads); break;
- case 2: gmshFunctionSpaces::find(MSH_TRI_6).df(uu,vv,grads); break;
- case 3: gmshFunctionSpaces::find(MSH_TRI_9).df(uu,vv,grads); break;
- case 4: gmshFunctionSpaces::find(MSH_TRI_12).df(uu,vv,grads); break;
- case 5: gmshFunctionSpaces::find(MSH_TRI_15I).df(uu,vv,grads); break;
+ case 1: gmshFunctionSpaces::find(MSH_TRI_3).df(uu, vv, grads); break;
+ case 2: gmshFunctionSpaces::find(MSH_TRI_6).df(uu, vv, grads); break;
+ case 3: gmshFunctionSpaces::find(MSH_TRI_9).df(uu, vv, grads); break;
+ case 4: gmshFunctionSpaces::find(MSH_TRI_12).df(uu, vv, grads); break;
+ case 5: gmshFunctionSpaces::find(MSH_TRI_15I).df(uu, vv, grads); break;
default: throw;
}
}
else{
switch(ord){
- case 1: gmshFunctionSpaces::find(MSH_TRI_3).df(uu,vv,grads); break;
- case 2: gmshFunctionSpaces::find(MSH_TRI_6).df(uu,vv,grads); break;
- case 3: gmshFunctionSpaces::find(MSH_TRI_10).df(uu,vv,grads); break;
- case 4: gmshFunctionSpaces::find(MSH_TRI_15).df(uu,vv,grads); break;
- case 5: gmshFunctionSpaces::find(MSH_TRI_21).df(uu,vv,grads); break;
+ case 1: gmshFunctionSpaces::find(MSH_TRI_3).df(uu, vv, grads); break;
+ case 2: gmshFunctionSpaces::find(MSH_TRI_6).df(uu, vv, grads); break;
+ case 3: gmshFunctionSpaces::find(MSH_TRI_10).df(uu, vv, grads); break;
+ case 4: gmshFunctionSpaces::find(MSH_TRI_15).df(uu, vv, grads); break;
+ case 5: gmshFunctionSpaces::find(MSH_TRI_21).df(uu, vv, grads); break;
default: throw;
}
}
- j[0][0] = 0 ; for(int i = 0; i < 3; i++) j[0][0] += grads [i][0] * _v[i] -> x();
- j[1][0] = 0 ; for(int i = 0; i < 3; i++) j[1][0] += grads [i][1] * _v[i] -> x();
- j[0][1] = 0 ; for(int i = 0; i < 3; i++) j[0][1] += grads [i][0] * _v[i] -> y();
- j[1][1] = 0 ; for(int i = 0; i < 3; i++) j[1][1] += grads [i][1] * _v[i] -> y();
- j[0][2] = 0 ; for(int i = 0; i < 3; i++) j[0][2] += grads [i][0] * _v[i] -> z();
- j[1][2] = 0 ; for(int i = 0; i < 3; i++) j[1][2] += grads [i][1] * _v[i] -> z();
-
- if (ord == 1)return;
-
- for(int i = 3; i < 3 * ord + nf; i++) j[0][0] += grads[i][0] * vs[i - 3] -> x();
- for(int i = 3; i < 3 * ord + nf; i++) j[1][0] += grads[i][1] * vs[i - 3] -> x();
- for(int i = 3; i < 3 * ord + nf; i++) j[0][1] += grads[i][0] * vs[i - 3] -> y();
- for(int i = 3; i < 3 * ord + nf; i++) j[1][1] += grads[i][1] * vs[i - 3] -> y();
- for(int i = 3; i < 3 * ord + nf; i++) j[0][2] += grads[i][0] * vs[i - 3] -> z();
- for(int i = 3; i < 3 * ord + nf; i++) j[1][2] += grads[i][1] * vs[i - 3] -> z();
+ j[0][0] = 0 ; for(int i = 0; i < 3; i++) j[0][0] += grads [i][0] * _v[i]->x();
+ j[1][0] = 0 ; for(int i = 0; i < 3; i++) j[1][0] += grads [i][1] * _v[i]->x();
+ j[0][1] = 0 ; for(int i = 0; i < 3; i++) j[0][1] += grads [i][0] * _v[i]->y();
+ j[1][1] = 0 ; for(int i = 0; i < 3; i++) j[1][1] += grads [i][1] * _v[i]->y();
+ j[0][2] = 0 ; for(int i = 0; i < 3; i++) j[0][2] += grads [i][0] * _v[i]->z();
+ j[1][2] = 0 ; for(int i = 0; i < 3; i++) j[1][2] += grads [i][1] * _v[i]->z();
+
+ if (ord == 1) return;
+
+ for(int i = 3; i < 3 * ord + nf; i++) j[0][0] += grads[i][0] * vs[i - 3]->x();
+ for(int i = 3; i < 3 * ord + nf; i++) j[1][0] += grads[i][1] * vs[i - 3]->x();
+ for(int i = 3; i < 3 * ord + nf; i++) j[0][1] += grads[i][0] * vs[i - 3]->y();
+ for(int i = 3; i < 3 * ord + nf; i++) j[1][1] += grads[i][1] * vs[i - 3]->y();
+ for(int i = 3; i < 3 * ord + nf; i++) j[0][2] += grads[i][0] * vs[i - 3]->z();
+ for(int i = 3; i < 3 * ord + nf; i++) j[1][2] += grads[i][1] * vs[i - 3]->z();
+#endif
}
void MTriangle::pnt(int ord, MVertex *vs[], double uu, double vv, SPoint3 &p)
{
+#if !defined(HAVE_GMSH_EMBEDDED)
double sf[256];
int nf = getNumFaceVertices();
if (!nf){
switch(ord){
- case 1: gmshFunctionSpaces::find(MSH_TRI_3).f(uu,vv,sf); break;
- case 2: gmshFunctionSpaces::find(MSH_TRI_6).f(uu,vv,sf); break;
- case 3: gmshFunctionSpaces::find(MSH_TRI_9).f(uu,vv,sf); break;
- case 4: gmshFunctionSpaces::find(MSH_TRI_12).f(uu,vv,sf); break;
- case 5: gmshFunctionSpaces::find(MSH_TRI_15I).f(uu,vv,sf); break;
+ case 1: gmshFunctionSpaces::find(MSH_TRI_3).f(uu, vv, sf); break;
+ case 2: gmshFunctionSpaces::find(MSH_TRI_6).f(uu, vv, sf); break;
+ case 3: gmshFunctionSpaces::find(MSH_TRI_9).f(uu, vv, sf); break;
+ case 4: gmshFunctionSpaces::find(MSH_TRI_12).f(uu, vv, sf); break;
+ case 5: gmshFunctionSpaces::find(MSH_TRI_15I).f(uu, vv, sf); break;
default: throw;
}
}
else{
switch(ord){
- case 1: gmshFunctionSpaces::find(MSH_TRI_3).f(uu,vv,sf); break;
- case 2: gmshFunctionSpaces::find(MSH_TRI_6).f(uu,vv,sf); break;
- case 3: gmshFunctionSpaces::find(MSH_TRI_10).f(uu,vv,sf); break;
- case 4: gmshFunctionSpaces::find(MSH_TRI_15).f(uu,vv,sf); break;
- case 5: gmshFunctionSpaces::find(MSH_TRI_21).f(uu,vv,sf); break;
+ case 1: gmshFunctionSpaces::find(MSH_TRI_3).f(uu, vv, sf); break;
+ case 2: gmshFunctionSpaces::find(MSH_TRI_6).f(uu, vv, sf); break;
+ case 3: gmshFunctionSpaces::find(MSH_TRI_10).f(uu, vv, sf); break;
+ case 4: gmshFunctionSpaces::find(MSH_TRI_15).f(uu, vv, sf); break;
+ case 5: gmshFunctionSpaces::find(MSH_TRI_21).f(uu, vv, sf); break;
default: throw;
}
}
- double x = 0 ; for(int i = 0; i < 3; i++) x += sf[i] * _v[i] -> x();
- double y = 0 ; for(int i = 0; i < 3; i++) y += sf[i] * _v[i] -> y();
- double z = 0 ; for(int i = 0; i < 3; i++) z += sf[i] * _v[i] -> z();
+ double x = 0 ; for(int i = 0; i < 3; i++) x += sf[i] * _v[i]->x();
+ double y = 0 ; for(int i = 0; i < 3; i++) y += sf[i] * _v[i]->y();
+ double z = 0 ; for(int i = 0; i < 3; i++) z += sf[i] * _v[i]->z();
- for(int i = 3; i < 3 * ord + nf; i++) x += sf[i] * vs[i - 3] -> x();
- for(int i = 3; i < 3 * ord + nf; i++) y += sf[i] * vs[i - 3] -> y();
- for(int i = 3; i < 3 * ord + nf; i++) z += sf[i] * vs[i - 3] -> z();
+ for(int i = 3; i < 3 * ord + nf; i++) x += sf[i] * vs[i - 3]->x();
+ for(int i = 3; i < 3 * ord + nf; i++) y += sf[i] * vs[i - 3]->y();
+ for(int i = 3; i < 3 * ord + nf; i++) z += sf[i] * vs[i - 3]->z();
p = SPoint3(x,y,z);
-
+#endif
}
-
-
void MTriangleN::jac(double uu, double vv , double j[2][3])
{
MTriangle::jac(_order, &(*(_vs.begin())), uu, vv, j);
@@ -754,68 +756,69 @@ void MTriangle::pnt(double uu, double vv, SPoint3 &p){
}
int MTriangle6::getNumEdgesRep(){ return 30; }
-void MTriangle6::getEdgeRep (int num, double *x, double *y, double *z, SVector3 *n){
+
+void MTriangle6::getEdgeRep (int num, double *x, double *y, double *z, SVector3 *n)
+{
if (num < 10){
- SPoint3 pnt1,pnt2;
- pnt ( (double)num/10. , 0. , pnt1);
- pnt ( (double)(num+1)/10. , 0. , pnt2);
- x[0] = pnt1.x();x[1] = pnt2.x();
- y[0] = pnt1.y();y[1] = pnt2.y();
- z[0] = pnt1.z();z[1] = pnt2.z();
+ SPoint3 pnt1, pnt2;
+ pnt((double)num / 10., 0., pnt1);
+ pnt((double)(num + 1) / 10., 0., pnt2);
+ x[0] = pnt1.x(); x[1] = pnt2.x();
+ y[0] = pnt1.y(); y[1] = pnt2.y();
+ z[0] = pnt1.z(); z[1] = pnt2.z();
return;
}
-
if (num < 20){
- SPoint3 pnt1,pnt2;
- num -=10;
- pnt ( 1.-(double)num/10. , (double)num/10. , pnt1);
- pnt ( 1.-(double)(num+1)/10. , (double)(num+1)/10. , pnt2);
- x[0] = pnt1.x();x[1] = pnt2.x();
- y[0] = pnt1.y();y[1] = pnt2.y();
- z[0] = pnt1.z();z[1] = pnt2.z();
+ SPoint3 pnt1, pnt2;
+ num -= 10;
+ pnt(1. - (double)num / 10., (double)num / 10., pnt1);
+ pnt(1. - (double)(num + 1) / 10., (double)(num + 1) / 10., pnt2);
+ x[0] = pnt1.x(); x[1] = pnt2.x();
+ y[0] = pnt1.y(); y[1] = pnt2.y();
+ z[0] = pnt1.z(); z[1] = pnt2.z();
return ;
}
{
- SPoint3 pnt1,pnt2;
+ SPoint3 pnt1, pnt2;
num -= 20;
- pnt ( 0,(double)num/10. , pnt1);
- pnt ( 0,(double)(num+1)/10., pnt2);
- x[0] = pnt1.x();x[1] = pnt2.x();
- y[0] = pnt1.y();y[1] = pnt2.y();
- z[0] = pnt1.z();z[1] = pnt2.z();
+ pnt(0, (double)num / 10., pnt1);
+ pnt(0, (double)(num + 1) / 10., pnt2);
+ x[0] = pnt1.x(); x[1] = pnt2.x();
+ y[0] = pnt1.y(); y[1] = pnt2.y();
+ z[0] = pnt1.z(); z[1] = pnt2.z();
}
}
int MTriangleN::getNumEdgesRep(){ return 120; }
-void MTriangleN::getEdgeRep (int num, double *x, double *y, double *z, SVector3 *n){
+
+void MTriangleN::getEdgeRep (int num, double *x, double *y, double *z, SVector3 *n)
+{
if (num < 40){
- SPoint3 pnt1,pnt2;
- pnt ( (double)num/40. , 0. , pnt1);
- pnt ( (double)(num+1)/40. , 0. , pnt2);
- x[0] = pnt1.x();x[1] = pnt2.x();
- y[0] = pnt1.y();y[1] = pnt2.y();
- z[0] = pnt1.z();z[1] = pnt2.z();
+ SPoint3 pnt1, pnt2;
+ pnt((double)num / 40., 0., pnt1);
+ pnt((double)(num + 1) / 40., 0., pnt2);
+ x[0] = pnt1.x(); x[1] = pnt2.x();
+ y[0] = pnt1.y(); y[1] = pnt2.y();
+ z[0] = pnt1.z(); z[1] = pnt2.z();
return;
}
-
if (num < 80){
- SPoint3 pnt1,pnt2;
- num -=40;
- pnt ( 1.-(double)num/40. , (double)num/40. , pnt1);
- pnt ( 1.-(double)(num+1)/40. , (double)(num+1)/40. , pnt2);
- x[0] = pnt1.x();x[1] = pnt2.x();
- y[0] = pnt1.y();y[1] = pnt2.y();
- z[0] = pnt1.z();z[1] = pnt2.z();
+ SPoint3 pnt1, pnt2;
+ num -= 40;
+ pnt(1. - (double)num / 40., (double)num / 40., pnt1);
+ pnt(1. - (double)(num + 1) / 40., (double)(num + 1) / 40., pnt2);
+ x[0] = pnt1.x(); x[1] = pnt2.x();
+ y[0] = pnt1.y(); y[1] = pnt2.y();
+ z[0] = pnt1.z(); z[1] = pnt2.z();
return ;
}
{
- SPoint3 pnt1,pnt2;
+ SPoint3 pnt1, pnt2;
num -= 80;
- pnt ( 0,(double)num/40. , pnt1);
- pnt ( 0,(double)(num+1)/40., pnt2);
- x[0] = pnt1.x();x[1] = pnt2.x();
- y[0] = pnt1.y();y[1] = pnt2.y();
- z[0] = pnt1.z();z[1] = pnt2.z();
+ pnt(0, (double)num / 40., pnt1);
+ pnt(0, (double)(num + 1) / 40., pnt2);
+ x[0] = pnt1.x(); x[1] = pnt2.x();
+ y[0] = pnt1.y(); y[1] = pnt2.y();
+ z[0] = pnt1.z(); z[1] = pnt2.z();
}
}
-
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 87adcd2fdf2d9da52ea61465d3c6d3e36a539322..597f9938308f5917149efef033637895592c785d 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -411,14 +411,14 @@ class MTriangle6 : public MTriangle {
virtual int getNumEdgeVertices(){ return 3; }
virtual int getNumEdgesRep();
virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n);
-/* { */
-/* static const int e[6][2] = { */
-/* {0, 3}, {3, 1}, */
-/* {1, 4}, {4, 2}, */
-/* {2, 5}, {5, 0} */
-/* }; */
-/* _getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n, 0); */
-/* } */
+ //{
+ // static const int e[6][2] = {
+ // {0, 3}, {3, 1},
+ // {1, 4}, {4, 2},
+ // {2, 5}, {5, 0}
+ // };
+ // _getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n, 0);
+ //}
virtual int getNumFacesRep(){ return 4; }
virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
{
diff --git a/Geo/Makefile b/Geo/Makefile
index db2abdb69dba87fe1275cc50f15b97a9264512ac..d1a24568686eb7890f2bc8c1960a52c1faf4f538 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.186 2008-02-21 09:49:22 geuzaine Exp $
+# $Id: Makefile,v 1.187 2008-02-21 12:11:12 geuzaine Exp $
#
# Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
#
@@ -82,16 +82,16 @@ GVertex.o: GVertex.cpp GVertex.h GEntity.h Range.h SPoint3.h \
GEdge.o: GEdge.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h SVector3.h GFace.h \
GEdgeLoop.h Pair.h GRegion.h MElement.h ../Common/GmshDefines.h \
- MVertex.h MEdge.h MFace.h ../Numeric/GaussLegendre1D.h \
- ../Common/Message.h
+ MVertex.h MEdge.h MFace.h ../Common/Message.h \
+ ../Numeric/GaussLegendre1D.h
GEdgeLoop.o: GEdgeLoop.cpp GEdgeLoop.h GEdge.h GEntity.h Range.h \
SPoint3.h SBoundingBox3d.h GVertex.h GPoint.h SPoint2.h SVector3.h \
../Common/Message.h
GFace.o: GFace.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h SVector3.h GFace.h \
GEdgeLoop.h Pair.h GRegion.h MElement.h ../Common/GmshDefines.h \
- MVertex.h MEdge.h MFace.h ../Numeric/GaussLegendre1D.h \
- ../Common/Message.h ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
+ MVertex.h MEdge.h MFace.h ../Common/Message.h ../Numeric/Numeric.h \
+ ../Numeric/NumericEmbedded.h ../Numeric/GaussLegendre1D.h \
../Common/VertexArray.h ../Geo/SVector3.h ../Common/Context.h \
../DataStr/List.h
GRegion.o: GRegion.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
@@ -246,6 +246,7 @@ MFace.o: MFace.cpp MFace.h MVertex.h SPoint3.h SVector3.h \
MElement.o: MElement.cpp MElement.h ../Common/GmshDefines.h MVertex.h \
SPoint3.h MEdge.h SVector3.h MFace.h GEntity.h Range.h SBoundingBox3d.h \
GFace.h GPoint.h GEdgeLoop.h GEdge.h GVertex.h SPoint2.h Pair.h \
- ../Numeric/FunctionSpace.h ../Common/GmshMatrix.h ../Numeric/Numeric.h \
- ../Numeric/NumericEmbedded.h ../Common/Message.h ../Common/Context.h \
- ../DataStr/List.h ../Mesh/qualityMeasures.h
+ ../Common/GmshMatrix.h ../Numeric/Numeric.h \
+ ../Numeric/NumericEmbedded.h ../Numeric/FunctionSpace.h \
+ ../Common/Message.h ../Common/Context.h ../DataStr/List.h \
+ ../Mesh/qualityMeasures.h
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 6817bdca9731a63c169459af0e974ef526ea828c..1d92fdffdbc5d65ac557fe36b81afd93ebad7c9c 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1,4 +1,4 @@
-// $Id: BDS.cpp,v 1.101 2008-02-21 09:45:15 remacle Exp $
+// $Id: BDS.cpp,v 1.102 2008-02-21 12:11:12 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -1248,8 +1248,8 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
double U = 0;
double V = 0;
double LC = 0;
- double oldU=p->u;
- double oldV=p->v;
+ double oldU = p->u;
+ double oldV = p->v;
std::list<BDS_Face*> ts;
p->getTriangles(ts);
@@ -1290,7 +1290,7 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
std::list<BDS_Face*>::iterator it = ts.begin();
std::list<BDS_Face*>::iterator ite = ts.end();
- double s1=0,s2=0;
+ double s1 = 0, s2 = 0;
double newWorst = 1.0;
double oldWorst = 1.0;
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 1e90129e24305e9fd9f67f9742161e1ef42c4313..8bb04cdcf9a5b93e189651c23351ce11090a2c8c 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -1,4 +1,4 @@
-// $Id: BackgroundMesh.cpp,v 1.37 2008-02-21 09:45:15 remacle Exp $
+// $Id: BackgroundMesh.cpp,v 1.38 2008-02-21 12:11:12 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -100,7 +100,6 @@ static double max_surf_curvature(const GVertex *gv)
return max_curvature;
}
-
// the mesh vertex is classified on a model vertex. we compute the
// maximum of the curvature of model faces surrounding this point if
// it is classified on a model edge, we do the same for all model
@@ -112,17 +111,15 @@ double LC_MVertex_CURV(GEntity *ge, double U, double V)
double Crv = 0;
switch(ge->dim()){
case 0:
- //Crv = max_edge_curvature ( (const GVertex *)ge);
- // Crv = std::max(max_surf_curvature ( (const GVertex *)ge),Crv);
- Crv = max_surf_curvature ( (const GVertex *)ge);
- // printf("point %d coucou %g\n",ge->tag(),Crv);
+ // Crv = max_edge_curvature((const GVertex *)ge);
+ // Crv = std::max(max_surf_curvature((const GVertex *)ge), Crv);
+ Crv = max_surf_curvature((const GVertex *)ge);
break;
case 1:
{
GEdge *ged = (GEdge *)ge;
- // Crv = ged->curvature(U);
- // printf("coucou %12.5E %d\n",Crv,CTX.mesh.min_circ_points);
- // Crv = std::max(Crv,max_surf_curvature(ged, U));
+ // Crv = ged->curvature(U);
+ // Crv = std::max(Crv, max_surf_curvature(ged, U));
Crv = max_surf_curvature(ged, U);
}
break;
@@ -137,13 +134,9 @@ double LC_MVertex_CURV(GEntity *ge, double U, double V)
double lc = Crv > 0 ? 2*M_PI / Crv / CTX.mesh.min_circ_points : MAX_LC;
// double lc_min = CTX.lc /300;
-
-
return lc;
-
}
-
// compute the mesh size at a given vertex due to prescribed sizes at
// mesh vertices
double LC_MVertex_PNTS(GEntity *ge, double U, double V)
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 163b139193dd8bd6b69f93679ddf46a5f29ed312..a6f04061479fa4efd48f810ab30180a20720f80a 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1,4 +1,4 @@
-// $Id: HighOrder.cpp,v 1.20 2008-02-21 09:45:15 remacle Exp $
+// $Id: HighOrder.cpp,v 1.21 2008-02-21 12:11:12 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -73,29 +73,26 @@ bool reparamOnFace(MVertex *v, GFace *gf, SPoint2 ¶m)
return true;
}
+// The aim here is to build a polynomial representation that consist
+// in polynomial segments of equal length
-/*
- The aim here is to build a polynomial representation
- that consist in polynomial segments of equal length
-*/
-
-static double mylength ( GEdge *ge , int i, double *u){
- // printf("from %22.15E to %22.15E (curve %d)\n",u[i],u[i+1],ge->tag());
- return ge->length ( u[i] , u[i+1], 6 );
+static double mylength(GEdge *ge, int i, double *u)
+{
+ return ge->length(u[i], u[i+1], 6);
}
-static void myresid ( int N , GEdge *ge , double *u, Double_Vector &r ){
+static void myresid(int N, GEdge *ge, double *u, Double_Vector &r)
+{
double L[100];
- for (int i=0;i<N-1;i++)L[i] = mylength(ge,i,u);
- for (int i=0;i<N-2;i++)r(i) = L[i+1]-L[i];
- // printf("%22.15E %22.15E %22.15E\n",L[0],L[1],r(0));
+ for (int i = 0; i < N - 1; i++) L[i] = mylength(ge, i, u);
+ for (int i = 0; i < N - 2; i++) r(i) = L[i + 1] - L[i];
}
-bool computeEquidistantParameters ( GEdge *ge , double u0, double uN , int N, double *u){
-#define NR_PRECISION 1.e-6
-#define NR_MAX_ITER 50
-
- const double eps = (uN-u0)*1.e-5;
+bool computeEquidistantParameters(GEdge *ge, double u0, double uN, int N, double *u)
+{
+ const double PRECISION = 1.e-6;
+ const int MAX_ITER = 50;
+ const double eps = (uN - u0) * 1.e-5;
// newton algorithm
// N is the total number of points (3 for quadratic, 4 for cubic ...)
@@ -105,78 +102,73 @@ bool computeEquidistantParameters ( GEdge *ge , double u0, double uN , int N, do
// initialize as equidistant in parameter space
u[0] = u0;
- double du = (uN-u0)/(N-1);
- for (int i=1 ; i < N ; i++){
- u[i] = u[i-1] + du;
+ double du = (uN - u0) / (N - 1);
+ for (int i = 1 ; i < N ; i++){
+ u[i] = u[i - 1] + du;
}
- // create the tangent matrix
- const int M = N-2;
- Double_Matrix J(M,M);
+ // create the tangent matrix
+ const int M = N - 2;
+ Double_Matrix J(M, M);
Double_Vector DU(M);
Double_Vector R(M);
Double_Vector Rp(M);
- int iter = 1 ;
+ int iter = 1;
double L[100];
- // printf("M = %d u = %12.5E %12.5E %12.5E\n",M,u[0],u[1],u[2]);
-
- while (iter < NR_MAX_ITER){
-
- iter++ ;
+ while (iter < MAX_ITER){
+ iter++;
+ myresid(N, ge, u, R);
- myresid ( N, ge, u , R);
-
- for (int i=0;i<M;i++){
- u[i+1] += eps;
- myresid ( N, ge, u, Rp);
- for (int j=0;j<M;j++){
- J(i,j) = ( Rp(j) - R(j)) / eps;
+ for (int i = 0; i < M; i++){
+ u[i + 1] += eps;
+ myresid(N, ge, u, Rp);
+ for (int j = 0; j < M; j++){
+ J(i, j) = (Rp(j) - R(j)) / eps;
}
u[i+1] -= eps;
}
-
- // printf("J(0,0) = %12.5E\n",J(0,0));
+
if (M == 1)
- DU(0) = R(0)/J(0,0);
+ DU(0) = R(0) / J(0, 0);
else
- J.lu_solve(R,DU);
+ J.lu_solve(R, DU);
- for (int i=0;i<M;i++){
+ for (int i = 0; i < M; i++){
u[i+1] -= DU(i);
}
- if (u[1] < u0)break;
- if (u[N-1] > uN)break;
+ if (u[1] < u0) break;
+ if (u[N - 1] > uN) break;
double newt_norm = DU.norm();
- // printf("iter %d u = %12.5E %12.5E %12.5E res %12.5E du %12.5E\n",iter,u[0],u[1],u[2],R(0),DU(0));
-
- if (newt_norm < NR_PRECISION)return true;
+ if (newt_norm < PRECISION) return true;
}
// FAILED, use equidistant in param space
- for (int i=1 ; i < N ; i++){
- u[i] = u[i-1] + du;
+ for (int i = 1; i < N; i++){
+ u[i] = u[i - 1] + du;
}
return false;
}
-static double mylength ( GFace *gf , int i, double *u, double *v){
- return gf->length ( SPoint2(u[i],v[i]) , SPoint2(u[i+1],v[i+1]), 4 );
+static double mylength(GFace *gf, int i, double *u, double *v)
+{
+ return gf->length(SPoint2(u[i], v[i]), SPoint2(u[i + 1], v[i + 1]), 4);
}
-static void myresid ( int N , GFace *gf , double *u, double *v, Double_Vector &r ){
+static void myresid(int N, GFace *gf, double *u, double *v, Double_Vector &r)
+{
double L[100];
- for (int i=0;i<N-1;i++)L[i] = mylength(gf,i,u,v);
- for (int i=0;i<N-2;i++)r(i) = L[i+1]-L[i];
- // printf("%22.15E %22.15E %22.15E\n",L[0],L[1],r(0));
+ for (int i = 0; i < N - 1; i++) L[i] = mylength(gf, i, u, v);
+ for (int i = 0; i < N - 2; i++) r(i) = L[i + 1] - L[i];
}
-bool computeEquidistantParameters ( GFace *gf , double u0, double uN , double v0, double vN, int N, double *u, double *v){
-#define NR_PRECISION 1.e-6
-#define NR_MAX_ITER 50
-
+bool computeEquidistantParameters(GFace *gf, double u0, double uN, double v0, double vN,
+ int N, double *u, double *v)
+{
+ const double PRECISION = 1.e-6;
+ const int MAX_ITER = 50;
const double eps = 1.e-4;
double t[100];
@@ -185,9 +177,9 @@ bool computeEquidistantParameters ( GFace *gf , double u0, double uN , double v0
u[0] = u0;
v[0] = v0;
t[0] = 0;
- for (int i=1 ; i < N ; i++){
- t[i] = (double)i/(N-1);
- SPoint2 p = gf->geodesic ( SPoint2(u0,v0), SPoint2(uN,vN), t[i] );
+ for (int i = 1; i < N ; i++){
+ t[i] = (double)i / (N - 1);
+ SPoint2 p = gf->geodesic(SPoint2(u0, v0), SPoint2(uN, vN), t[i]);
u[i] = p.x();
v[i] = p.y();
}
@@ -196,62 +188,54 @@ bool computeEquidistantParameters ( GFace *gf , double u0, double uN , double v0
t[N] = 1.0;
// create the tangent matrix
-
- const int M = N-2;
- Double_Matrix J(M,M);
+ const int M = N - 2;
+ Double_Matrix J(M, M);
Double_Vector DU(M);
Double_Vector R(M);
Double_Vector Rp(M);
- int iter = 1 ;
+ int iter = 1 ;
double L[100];
- // printf("M = %d u = %12.5E %12.5E %12.5E\n",M,u[0],u[1],u[2]);
-
- while (iter < NR_MAX_ITER){
-
+ while (iter < MAX_ITER){
iter++ ;
-
- myresid ( N, gf, u , v, R);
-
- for (int i=0;i<M;i++){
- t[i+1] += eps;
- double tempu = u[i+1];
- double tempv = v[i+1];
- SPoint2 p = gf->geodesic ( SPoint2(u0,v0), SPoint2(uN,vN), t[i+1] );
- u[i+1] = p.x();
- v[i+1] = p.y();
- myresid ( N, gf, u, v, Rp);
- for (int j=0;j<M;j++){
- J(i,j) = ( Rp(j) - R(j)) / eps;
+ myresid(N, gf, u, v, R);
+
+ for (int i = 0; i < M; i++){
+ t[i + 1] += eps;
+ double tempu = u[i + 1];
+ double tempv = v[i + 1];
+ SPoint2 p = gf->geodesic(SPoint2(u0, v0), SPoint2(uN, vN), t[i + 1]);
+ u[i + 1] = p.x();
+ v[i + 1] = p.y();
+ myresid(N, gf, u, v, Rp);
+ for (int j = 0; j < M; j++){
+ J(i, j) = (Rp(j) - R(j)) / eps;
}
- t[i+1] -= eps;
- u[i+1] = tempu;
- v[i+1] = tempv;
+ t[i + 1] -= eps;
+ u[i + 1] = tempu;
+ v[i + 1] = tempv;
}
-
- // printf("J(0,0) = %12.5E\n",J(0,0));
+
if (M == 1)
- DU(0) = R(0)/J(0,0);
+ DU(0) = R(0) / J(0, 0);
else
- J.lu_solve(R,DU);
+ J.lu_solve(R, DU);
- for (int i=0;i<M;i++){
- t[i+1] -= DU(i);
- SPoint2 p = gf->geodesic ( SPoint2(u0,v0), SPoint2(uN,vN), t[i+1] );
- u[i+1] = p.x();
- v[i+1] = p.y();
+ for (int i = 0; i < M; i++){
+ t[i + 1] -= DU(i);
+ SPoint2 p = gf->geodesic(SPoint2(u0, v0), SPoint2(uN, vN), t[i + 1]);
+ u[i + 1] = p.x();
+ v[i + 1] = p.y();
}
double newt_norm = DU.norm();
- // printf("iter %d u = %12.5E %12.5E %12.5E res %12.5E du %12.5E\n",iter,u[0],u[1],u[2],R(0),DU(0));
-
- if (newt_norm < NR_PRECISION)return true;
+ if (newt_norm < PRECISION) return true;
}
// FAILED, use equidistant in param space
- for (int i=1 ; i < N ; i++){
- t[i] = (double)i/(N-1);
- SPoint2 p = gf->geodesic ( SPoint2(u0,v0), SPoint2(uN,vN), t[i] );
+ for (int i = 1; i < N; i++){
+ t[i] = (double)i / (N - 1);
+ SPoint2 p = gf->geodesic(SPoint2(u0, v0), SPoint2(uN, vN), t[i]);
u[i] = p.x();
v[i] = p.y();
}
@@ -294,13 +278,11 @@ void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
reparamOK &= reparamOnEdge(v0, ge, u0);
reparamOK &= reparamOnEdge(v1, ge, u1);
}
- std::vector<MVertex*> temp;
-
double US[100];
if(reparamOK && !linear && ge->geomType() != GEntity::DiscreteCurve){
- computeEquidistantParameters ( ge , u0, u1 , nPts+2, US);
+ computeEquidistantParameters(ge, u0, u1, nPts + 2, US);
}
-
+ std::vector<MVertex*> temp;
for(int j = 0; j < nPts; j++){
const double t = (double)(j + 1)/(nPts + 1);
double uc = (1. - t) * u0 + t * u1;
@@ -349,13 +331,11 @@ void getEdgeVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &ve,
reparamOK &= reparamOnFace(v0, gf, p0);
reparamOK &= reparamOnFace(v1, gf, p1);
}
- std::vector<MVertex*> temp;
-
- double US[100],VS[100];
+ double US[100], VS[100];
if(reparamOK && !linear && gf->geomType() != GEntity::DiscreteCurve){
- computeEquidistantParameters (gf , p0[0], p1[0], p0[1], p1[1], nPts+2, US,VS);
+ computeEquidistantParameters(gf, p0[0], p1[0], p0[1], p1[1], nPts + 2, US, VS);
}
-
+ std::vector<MVertex*> temp;
for(int j = 0; j < nPts; j++){
const double t = (double)(j + 1) / (nPts + 1);
MVertex *v;
@@ -418,15 +398,15 @@ void getFaceVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &vf,
int start;
switch (nPts){
case 2 :
- points= gmshFunctionSpaces::find(MSH_TRI_10).points;
+ points = gmshFunctionSpaces::find(MSH_TRI_10).points;
start = 9;
break;
case 3 :
- points= gmshFunctionSpaces::find(MSH_TRI_15).points;
+ points = gmshFunctionSpaces::find(MSH_TRI_15).points;
start = 12;
break;
case 4 :
- points= gmshFunctionSpaces::find(MSH_TRI_21).points;
+ points = gmshFunctionSpaces::find(MSH_TRI_21).points;
start = 15;
break;
default :
@@ -453,14 +433,10 @@ void getFaceVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &vf,
reparamOK &= reparamOnFace(ele->getVertex(3), gf, p3);
}
if(face.getNumVertices() == 3){ // triangles
- // for(int j = 0; j < nPts; j++){
- // for(int k = 0 ; k < nPts - j - 1; k++){
for(int k = start ; k < points.size1() ; k++){
MVertex *v;
- // double t1 = (double)(j + 1) / (nPts + 1);
- // double t2 = (double)(k + 1) / (nPts + 1);
- const double t1 = points(k,0);
- const double t2 = points(k,1);
+ const double t1 = points(k, 0);
+ const double t2 = points(k, 1);
if(!reparamOK || linear || gf->geomType() == GEntity::DiscreteSurface){
SPoint3 pc = face.interpolate(t1, t2);
v = new MVertex(pc.x(), pc.y(), pc.z(), gf);
@@ -797,40 +773,26 @@ bool straightLine(std::vector<MVertex*> &l, MVertex *n1, MVertex *n2)
return true;
}
-FILE *MYFILE = 0;
-
void getMinMaxJac (MTriangle *t, double &minJ, double &maxJ)
{
double mat[2][3];
int n = 3;
-
- t->jac(1,0,0,0,mat);
- double v1[3] = {mat[0][0],mat[0][1],mat[0][2]};
- double v2[3] = {mat[1][0],mat[1][1],mat[1][2]};
- double normal1[3],normal[3];
- prodve(v1,v2,normal1);
- // norme(normal1);
- double nn = sqrt(DSQR(normal1[0]) + DSQR(normal1[1]) + DSQR(normal1[2])) ;
-
- // double sign = mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0];
- // double invsign = 1./(sign);
-
+ t->jac(1, 0, 0, 0, mat);
+ double v1[3] = {mat[0][0], mat[0][1], mat[0][2]};
+ double v2[3] = {mat[1][0], mat[1][1], mat[1][2]};
+ double normal1[3], normal[3];
+ prodve(v1, v2, normal1);
+ double nn = sqrt(DSQR(normal1[0]) + DSQR(normal1[1]) + DSQR(normal1[2]));
for(int i = 0; i < n; i++){
for(int k = 0; k < n - i; k++){
t->jac((double)i / (n - 1), (double)k / (n - 1), mat);
- // printf("%g %g\n",(double)i / (n - 1), (double)k / (n - 1));
- // SPoint3 pt;
- // t->pnt((double)i / (n - 1), (double)k / (n - 1),pt);
- double v1b[3] = {mat[0][0],mat[0][1],mat[0][2]};
- double v2b[3] = {mat[1][0],mat[1][1],mat[1][2]};
- prodve(v1b,v2b,normal);
- double sign; prosca(normal1,normal,&sign);
- double det = norm3(normal) * (sign>0?1.:-1.) / nn;
- //double det = invsign*(mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0]) * (sign>0?1.:-1.);
- // double det2 = 1./det1;
- // double det = std::max(det1,det2);
- // if(MYFILE)fprintf(MYFILE,"SP(%g,%g,%g){%g};\n",pt.x(),pt.y(),pt.z(),det);
- minJ = std::min(1./det,std::min(det, minJ));
+ double v1b[3] = {mat[0][0], mat[0][1], mat[0][2]};
+ double v2b[3] = {mat[1][0], mat[1][1], mat[1][2]};
+ prodve(v1b, v2b, normal);
+ double sign;
+ prosca(normal1, normal, &sign);
+ double det = norm3(normal) * (sign > 0 ? 1. : -1.) / nn;
+ minJ = std::min(1. / det, std::min(det, minJ));
maxJ = std::max(det, maxJ);
}
}
@@ -839,18 +801,19 @@ void getMinMaxJac (MTriangle *t, double &minJ, double &maxJ)
struct smoothVertexDataHO{
MVertex *v;
GFace *gf;
- std::vector < MTriangle * >ts;
+ std::vector<MTriangle*> ts;
};
struct smoothVertexDataHON{
std::vector<MVertex*> v;
GFace *gf;
- std::vector < MTriangle * >ts;
+ std::vector<MTriangle*> ts;
};
-double smoothing_objective_function_HighOrder (double U, double V, MVertex *v, std::vector < MTriangle * >&ts ,GFace *gf){
-
- GPoint gp = gf->point(U,V);
+double smoothing_objective_function_HighOrder(double U, double V, MVertex *v,
+ std::vector<MTriangle*> &ts, GFace *gf)
+{
+ GPoint gp = gf->point(U, V);
const double oldX = v->x();
const double oldY = v->y();
const double oldZ = v->z();
@@ -872,28 +835,32 @@ double smoothing_objective_function_HighOrder (double U, double V, MVertex *v, s
}
void deriv_smoothing_objective_function_HighOrder(double U, double V,
- double &F, double &dFdU, double &dFdV,void *data){
+ double &F, double &dFdU,
+ double &dFdV, void *data)
+{
smoothVertexDataHO *svd = (smoothVertexDataHO*)data;
MVertex *v = svd->v;
const double LARGE = -1.e5;
const double SMALL = 1./LARGE;
- F = smoothing_objective_function_HighOrder(U,V,v,svd->ts,svd->gf);
- double F_U = smoothing_objective_function_HighOrder(U+SMALL,V,v,svd->ts,svd->gf);
- double F_V = smoothing_objective_function_HighOrder(U,V+SMALL,v,svd->ts,svd->gf);
- dFdU = (F_U-F)*LARGE;
- dFdV = (F_V-F)*LARGE;
+ F = smoothing_objective_function_HighOrder(U, V, v, svd->ts, svd->gf);
+ double F_U = smoothing_objective_function_HighOrder(U + SMALL, V, v, svd->ts, svd->gf);
+ double F_V = smoothing_objective_function_HighOrder(U, V + SMALL, v, svd->ts, svd->gf);
+ dFdU = (F_U - F) * LARGE;
+ dFdV = (F_V - F) * LARGE;
}
-double smooth_obj_HighOrder(double U, double V, void *data){
+double smooth_obj_HighOrder(double U, double V, void *data)
+{
smoothVertexDataHO *svd = (smoothVertexDataHO*)data;
- return smoothing_objective_function_HighOrder(U,V,svd->v,svd->ts,svd->gf);
+ return smoothing_objective_function_HighOrder(U, V, svd->v, svd->ts, svd->gf);
}
-double smooth_obj_HighOrderN(double *uv, void *data){
+double smooth_obj_HighOrderN(double *uv, void *data)
+{
smoothVertexDataHON *svd = (smoothVertexDataHON*)data;
double oldX[10],oldY[10],oldZ[10];
for (unsigned int i = 0; i < svd->v.size(); i++){
- GPoint gp = svd->gf->point(uv[2*i],uv[2*i+1]);
+ GPoint gp = svd->gf->point(uv[2 * i], uv[2 * i + 1]);
oldX[i] = svd->v[i]->x();
oldY[i] = svd->v[i]->y();
oldZ[i] = svd->v[i]->z();
@@ -914,42 +881,43 @@ double smooth_obj_HighOrderN(double *uv, void *data){
return -minJ;
}
-
-void deriv_smoothing_objective_function_HighOrderN(double *uv, double *dF, double &F,void *data){
+void deriv_smoothing_objective_function_HighOrderN(double *uv, double *dF,
+ double &F, void *data)
+{
const double LARGE = -1.e2;
- const double SMALL = 1./LARGE;
+ const double SMALL = 1. / LARGE;
smoothVertexDataHON *svd = (smoothVertexDataHON*)data;
- F = smooth_obj_HighOrderN(uv,data);
+ F = smooth_obj_HighOrderN(uv, data);
for (unsigned int i = 0; i < svd->v.size(); i++){
uv[i] += SMALL;
- dF[i] = (smooth_obj_HighOrderN(uv,data) - F)*LARGE;
+ dF[i] = (smooth_obj_HighOrderN(uv, data) - F) * LARGE;
uv[i] -= SMALL;
}
}
-void optimizeNodeLocations ( GFace *gf, smoothVertexDataHON &vdN , double eps = .2){
-
- if(!vdN.v.size())return;
+void optimizeNodeLocations ( GFace *gf, smoothVertexDataHON &vdN , double eps = .2)
+{
+ if(!vdN.v.size()) return;
double uv[20];
for (unsigned int i = 0; i < vdN.v.size(); i++){
- if (!vdN.v[i]->getParameter (0,uv[2*i]))throw;
- if (!vdN.v[i]->getParameter (1,uv[2*i+1]))throw;
+ if (!vdN.v[i]->getParameter(0, uv[2 * i])) throw;
+ if (!vdN.v[i]->getParameter(1, uv[2 * i + 1])) throw;
}
double F = -smooth_obj_HighOrderN(uv, &vdN);
if (F < eps){
double val;
- minimize_N ( 2*vdN.v.size(),
- smooth_obj_HighOrderN,
- deriv_smoothing_objective_function_HighOrderN,
- &vdN, 1, uv,val);
+ minimize_N(2 * vdN.v.size(),
+ smooth_obj_HighOrderN,
+ deriv_smoothing_objective_function_HighOrderN,
+ &vdN, 1, uv,val);
double Fafter = -smooth_obj_HighOrderN(uv, &vdN);
- printf("%12.5E %12.5E\n",F,Fafter);
+ printf("%12.5E %12.5E\n", F, Fafter);
if (F < Fafter){
for (unsigned int i = 0; i < vdN.v.size(); i++){
- vdN.v[i]->setParameter ( 0,uv[2*i]);
- vdN.v[i]->setParameter ( 1,uv[2*i+1]);
- GPoint gp = gf->point(uv[2*i],uv[2*i+1]);
+ vdN.v[i]->setParameter(0, uv[2 * i]);
+ vdN.v[i]->setParameter(1, uv[2 * i + 1]);
+ GPoint gp = gf->point(uv[2 * i], uv[2 * i + 1]);
vdN.v[i]->x() = gp.x();
vdN.v[i]->y() = gp.y();
vdN.v[i]->z() = gp.z();
@@ -958,8 +926,8 @@ void optimizeNodeLocations ( GFace *gf, smoothVertexDataHON &vdN , double eps =
}
}
-void optimizeHighOrderMeshInternalNodes(GFace *gf){
-
+void optimizeHighOrderMeshInternalNodes(GFace *gf)
+{
for(unsigned int i = 0; i < gf->triangles.size(); i++){
MTriangle *t = gf->triangles[i];
smoothVertexDataHON vdN;
@@ -968,17 +936,14 @@ void optimizeHighOrderMeshInternalNodes(GFace *gf){
vdN.v.push_back(t->getVertex(j));
vdN.gf = gf;
vdN.ts.push_back(t);
- optimizeNodeLocations ( gf, vdN , .9);
+ optimizeNodeLocations(gf, vdN, .9);
}
-
}
-
-
-bool optimizeHighOrderMesh(GFace *gf, edgeContainer &edgeVertices){
-
+bool optimizeHighOrderMesh(GFace *gf, edgeContainer &edgeVertices)
+{
v2t_cont adjv;
- buildVertexToTriangle ( gf->triangles , adjv );
+ buildVertexToTriangle(gf->triangles, adjv);
typedef std::map<std::pair<MVertex*, MVertex*>, std::vector<MElement*> > edge2tris;
edge2tris e2t;
@@ -1012,7 +977,8 @@ bool optimizeHighOrderMesh(GFace *gf, edgeContainer &edgeVertices){
double val;
double F = -smooth_obj_HighOrder(initu,initv, &vd);
if (F < .2){
- minimize_2 ( smooth_obj_HighOrder, deriv_smoothing_objective_function_HighOrder, &vd, 1, initu,initv,val);
+ minimize_2 ( smooth_obj_HighOrder,
+ deriv_smoothing_objective_function_HighOrder, &vd, 1, initu,initv,val);
double Fafter = -smooth_obj_HighOrder(initu,initv, &vd);
if (F < Fafter){
success = true;
@@ -1051,7 +1017,7 @@ bool optimizeHighOrderMesh(GFace *gf, edgeContainer &edgeVertices){
vdN.ts.clear();
vdN.ts.push_back(t1);
vdN.ts.push_back(t2);
- optimizeNodeLocations ( gf, vdN );
+ optimizeNodeLocations(gf, vdN);
}
}
}
@@ -1088,25 +1054,25 @@ bool smoothInternalEdges(GFace *gf, edgeContainer &edgeVertices)
MVertex *n1 = t1->getOtherVertex(n2, n4);
MVertex *n3 = t2->getOtherVertex(n2, n4);
if(n1 < n2)
- e1 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n1, n2)];
+ e1 = edgeVertices[std::make_pair<MVertex*, MVertex*>(n1, n2)];
else
- e1 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n2, n1)];
+ e1 = edgeVertices[std::make_pair<MVertex*, MVertex*>(n2, n1)];
if(n2 < n3)
- e2 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n2, n3)];
+ e2 = edgeVertices[std::make_pair<MVertex*, MVertex*>(n2, n3)];
else
- e2 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n3, n2)];
+ e2 = edgeVertices[std::make_pair<MVertex*, MVertex*>(n3, n2)];
if(n3 < n4)
- e3 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n3, n4)];
+ e3 = edgeVertices[std::make_pair<MVertex*, MVertex*>(n3, n4)];
else
- e3 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n4, n3)];
+ e3 = edgeVertices[std::make_pair<MVertex*, MVertex*>(n4, n3)];
if(n4 < n1)
- e4 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n4, n1)];
+ e4 = edgeVertices[std::make_pair<MVertex*, MVertex*>(n4, n1)];
else
- e4 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n1, n4)];
+ e4 = edgeVertices[std::make_pair<MVertex*, MVertex*>(n1, n4)];
if(n2 < n4)
- e = edgeVertices[std::make_pair<MVertex*, MVertex*> (n2, n4)];
+ e = edgeVertices[std::make_pair<MVertex*, MVertex*>(n2, n4)];
else
- e = edgeVertices[std::make_pair<MVertex*, MVertex*> (n4, n2)];
+ e = edgeVertices[std::make_pair<MVertex*, MVertex*>(n4, n2)];
if((!straightLine(e1, n1, n2) || !straightLine(e2, n2, n3) ||
!straightLine(e3, n3, n4) || !straightLine(e4, n4, n1))){
@@ -1168,8 +1134,6 @@ bool smoothInternalEdges(GFace *gf, edgeContainer &edgeVertices)
getMinMaxJac (t1, minJloc, maxJloc);
getMinMaxJac (t2, minJloc, maxJloc);
- // printf("relax %d minJ %12.5E minjLoc %12.5E\n",k,minJ,minJloc);
-
if (minJloc > minJ){
kopt = k;
minJ = minJloc;
@@ -1213,23 +1177,21 @@ void checkHighOrderTriangles(GModel *m)
MTriangle *t = (*it)->triangles[i];
if(t->getPolynomialOrder() > 1 && t->getPolynomialOrder() < 6){
getMinMaxJac (t, minJloc, maxJloc);
- // printf("%p is %12.5E\n",t,minJloc);
minJ = std::min(minJ, minJloc);
maxJ = std::max(maxJ, maxJloc);
-// if(minJloc * maxJloc < 0)
-// Msg(WARNING, "Triangle %d (%d %d %d) has negative Jacobian (on gFace %d)", t->getNum(),
-// t->getVertex(0)->getNum(), t->getVertex(1)->getNum(), t->getVertex(2)->getNum(),(*it)->tag());
}
}
minJGlob = std::min(minJGlob,minJ);
maxJGlob = std::max(maxJGlob,maxJ);
}
- if (minJGlob >= 0)Msg(INFO, "Jacobian Range (%12.5E,%12.5E)",minJGlob, maxJGlob);
- else Msg(WARNING, "Jacobian Range (%12.5E,%12.5E)",minJGlob, maxJGlob);
+ if (minJGlob >= 0) Msg(INFO, "Jacobian Range (%12.5E,%12.5E)", minJGlob, maxJGlob);
+ else Msg(WARNING, "Jacobian Range (%12.5E,%12.5E)", minJGlob, maxJGlob);
}
void printJacobians(GModel *m, const char *nm)
{
+ return;
+
const int n = 22;
double D[n][n];
double X[n][n];
@@ -1279,9 +1241,6 @@ void printJacobians(GModel *m, const char *nm)
}
}
}
-
-
-
}
fprintf(f,"};\n");
fclose(f);
@@ -1322,7 +1281,7 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
setHighOrder(*it, edgeVertices, faceVertices, linear, incomplete, nPts);
- printJacobians(m,"detjIni.pos");
+ printJacobians(m, "detjIni.pos");
if(CTX.mesh.smooth_internal_edges){
checkHighOrderTriangles(m);
@@ -1341,13 +1300,10 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
}
}
}
- printJacobians(m,"detjOpt.pos");
- // MYFILE = fopen("jacs.pos","w");
- // fprintf(MYFILE,"View \"\"{\n");
+
+ printJacobians(m, "detjOpt.pos");
checkHighOrderTriangles(m);
- // fprintf(MYFILE,"};\n");
- // fclose(MYFILE);
- MYFILE=0;
+
double t2 = Cpu();
Msg(INFO, "Mesh second order complete (%g s)", t2 - t1);
Msg(STATUS1, "Mesh");
diff --git a/Mesh/HighOrder.h b/Mesh/HighOrder.h
index 0c067658777ec63e33a31f8c6a170fd472fb88ed..44b68e1101513c027f6768e4572c60fae4024962 100644
--- a/Mesh/HighOrder.h
+++ b/Mesh/HighOrder.h
@@ -24,6 +24,6 @@
void SetOrder1(GModel *m);
void SetOrderN(GModel *m, int order, bool linear=true, bool incomplete=false);
-void checkHighOrderTriangles ( GModel *m );
+void checkHighOrderTriangles(GModel *m);
#endif
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index 708b5fc2bb0af1f30af54f4b6934396792917cee..3feb8b99def5c25a1e126dc87a7c90fb7344977e 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceBDS.cpp,v 1.6 2008-02-21 09:45:15 remacle Exp $
+// $Id: meshGFaceBDS.cpp,v 1.7 2008-02-21 12:11:12 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -41,85 +41,86 @@ extern Context_T CTX;
inline double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2)
{
- const double dx = p1->X-p2->X;
- const double dy = p1->Y-p2->Y;
- const double dz = p1->Z-p2->Z;
- const double l = sqrt(dx*dx+dy*dy+dz*dz);
+ const double dx = p1->X - p2->X;
+ const double dy = p1->Y - p2->Y;
+ const double dz = p1->Z - p2->Z;
+ const double l = sqrt(dx * dx + dy * dy + dz * dz);
return l;
}
inline double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2, GFace *f,
double SCALINGU, double SCALINGV)
{
-
- GPoint GP = f->point (SPoint2(0.5*(p1->u+p2->u)*SCALINGU,0.5*(p1->v+p2->v)*SCALINGV));
-
- const double dx1 = p1->X-GP.x();
- const double dy1 = p1->Y-GP.y();
- const double dz1 = p1->Z-GP.z();
- const double l1 = sqrt(dx1*dx1+dy1*dy1+dz1*dz1);
- const double dx2 = p2->X-GP.x();
- const double dy2 = p2->Y-GP.y();
- const double dz2 = p2->Z-GP.z();
- const double l2 = sqrt(dx2*dx2+dy2*dy2+dz2*dz2);
- // printf("%g %g\n",l1,l2);
- return l1+l2;
+ GPoint GP = f->point(SPoint2(0.5 * (p1->u + p2->u) * SCALINGU,
+ 0.5 * (p1->v + p2->v) * SCALINGV));
+
+ const double dx1 = p1->X - GP.x();
+ const double dy1 = p1->Y - GP.y();
+ const double dz1 = p1->Z - GP.z();
+ const double l1 = sqrt(dx1 * dx1 + dy1 * dy1 + dz1 * dz1);
+ const double dx2 = p2->X - GP.x();
+ const double dy2 = p2->Y - GP.y();
+ const double dz2 = p2->Z - GP.z();
+ const double l2 = sqrt(dx2 * dx2 + dy2 * dy2 + dz2 * dz2);
+ return l1 + l2;
}
inline double computeEdgeMiddleCoord(BDS_Point *p1, BDS_Point *p2, GFace *f,
double SCALINGU, double SCALINGV)
{
-
if (f->geomType() == GEntity::Plane)
return 0.5;
- GPoint GP = f->point (SPoint2(0.5*(p1->u+p2->u)*SCALINGU,0.5*(p1->v+p2->v)*SCALINGV));
+ GPoint GP = f->point (SPoint2(0.5 * (p1->u + p2->u) * SCALINGU,
+ 0.5 * (p1->v + p2->v) * SCALINGV));
- const double dx1 = p1->X-GP.x();
- const double dy1 = p1->Y-GP.y();
- const double dz1 = p1->Z-GP.z();
- const double l1 = sqrt(dx1*dx1+dy1*dy1+dz1*dz1);
- const double dx2 = p2->X-GP.x();
- const double dy2 = p2->Y-GP.y();
- const double dz2 = p2->Z-GP.z();
- const double l2 = sqrt(dx2*dx2+dy2*dy2+dz2*dz2);
+ const double dx1 = p1->X - GP.x();
+ const double dy1 = p1->Y - GP.y();
+ const double dz1 = p1->Z - GP.z();
+ const double l1 = sqrt(dx1 * dx1 + dy1 * dy1 + dz1 * dz1);
+ const double dx2 = p2->X - GP.x();
+ const double dy2 = p2->Y - GP.y();
+ const double dz2 = p2->Z - GP.z();
+ const double l2 = sqrt(dx2 * dx2 + dy2 * dy2 + dz2 * dz2);
if (l1 > l2)
- return 0.25 * (l1+l2) / l1;
+ return 0.25 * (l1 + l2) / l1;
else
- return 0.25 * (3*l2-l1)/l2;
+ return 0.25 * (3 * l2 - l1) / l2;
}
-inline double computeEdgeLinearLength(BDS_Edge *e, GFace *f, double SCALINGU, double SCALINGV)
+inline double computeEdgeLinearLength(BDS_Edge *e, GFace *f,
+ double SCALINGU, double SCALINGV)
{
if (f->geomType() == GEntity::Plane)
return e->length();
else
- return computeEdgeLinearLength(e->p1, e->p2,f,SCALINGU,SCALINGV);
+ return computeEdgeLinearLength(e->p1, e->p2, f, SCALINGU, SCALINGV);
}
-double NewGetLc(BDS_Point *p){
- return Extend1dMeshIn2dSurfaces () ? std::min(p->lc(),p->lcBGM()) : p->lcBGM();
+double NewGetLc(BDS_Point *p)
+{
+ return Extend1dMeshIn2dSurfaces() ? std::min(p->lc(), p->lcBGM()) : p->lcBGM();
}
double NewGetLc(BDS_Edge *e, GFace *f, double SCALINGU, double SCALINGV)
{
- double linearLength = computeEdgeLinearLength(e,f,SCALINGU,SCALINGV);
+ double linearLength = computeEdgeLinearLength(e, f, SCALINGU, SCALINGV);
double l1 = NewGetLc(e->p1);
double l2 = NewGetLc(e->p2);
- return 2*linearLength / (l1 + l2);
+ return 2 * linearLength / (l1 + l2);
}
double NewGetLc(BDS_Point *p1,BDS_Point *p2, GFace *f, double su, double sv)
{
- double linearLength = computeEdgeLinearLength(p1,p2,f,su,sv);
+ double linearLength = computeEdgeLinearLength(p1, p2, f, su, sv);
double l1 = NewGetLc(p1);
double l2 = NewGetLc(p2);
- return 2*linearLength / (l1 + l2);
+ return 2 * linearLength / (l1 + l2);
}
-/*Size field accuracy*/
-void computeMeshSizeFieldAccuracy (GFace *gf, BDS_Mesh &m, double &avg, double &max_e, double &min_e, int &nE, int &GS)
+void computeMeshSizeFieldAccuracy(GFace *gf, BDS_Mesh &m, double &avg,
+ double &max_e, double &min_e, int &nE, int &GS)
{
std::list<BDS_Edge*>::iterator it = m.edges.begin();
avg=0;
@@ -131,18 +132,19 @@ void computeMeshSizeFieldAccuracy (GFace *gf, BDS_Mesh &m, double &avg, double &
double sqr2 = sqrt(2.);
while (it!= m.edges.end()){
if (!(*it)->deleted){
- double lone = NewGetLc ( *it,gf,m.scalingU,m.scalingV);
- if (lone > oneoversqr2 && lone < sqr2 ) GS++;
- avg+=lone>1 ? (1./lone) - 1. : lone - 1.;
- max_e = std::max(max_e,lone);
- min_e = std::min(min_e,lone);
+ double lone = NewGetLc(*it, gf, m.scalingU, m.scalingV);
+ if (lone > oneoversqr2 && lone < sqr2) GS++;
+ avg += lone >1 ? (1. / lone) - 1. : lone - 1.;
+ max_e = std::max(max_e, lone);
+ min_e = std::min(min_e, lone);
nE++;
}
++it;
}
}
-/*Element shapes*/
-void computeElementShapes (GFace *gf, BDS_Mesh &m, double &worst, double &avg, double &best, int &nT, int &nbGQ)
+
+void computeElementShapes(GFace *gf, BDS_Mesh &m, double &worst, double &avg,
+ double &best, int &nT, int &nbGQ)
{
std::list<BDS_Face*>::iterator it = m.triangles.begin();
worst = 1.e22;
@@ -152,11 +154,11 @@ void computeElementShapes (GFace *gf, BDS_Mesh &m, double &worst, double &avg, d
nbGQ = 0;
while (it!= m.triangles.end()){
if (!(*it)->deleted){
- double q = qmTriangle(*it,QMTRI_RHO);
- if (q>.9) nbGQ++;
- avg+=q;
- worst = std::min(worst,q);
- best = std::max(best,q);
+ double q = qmTriangle(*it, QMTRI_RHO);
+ if (q > .9) nbGQ++;
+ avg += q;
+ worst = std::min(worst, q);
+ best = std::max(best, q);
nT++;
}
++it;
@@ -164,9 +166,7 @@ void computeElementShapes (GFace *gf, BDS_Mesh &m, double &worst, double &avg, d
avg /= nT;
}
-/*
- SWAP TESTS i.e. tell if swap should be done
-*/
+// SWAP TESTS i.e. tell if swap should be done
bool edgeSwapTestAngle(BDS_Edge *e, double min_cos)
{
@@ -178,9 +178,9 @@ bool edgeSwapTestAngle(BDS_Edge *e, double min_cos)
f2->getNodes(n2);
double norm1[3];
double norm2[3];
- normal_triangle(n1[0],n1[1],n1[2],norm1);
- normal_triangle(n2[0],n2[1],n2[2],norm2);
- double cosa;prosca (norm1,norm2,&cosa);
+ normal_triangle(n1[0], n1[1], n1[2], norm1);
+ normal_triangle(n2[0], n2[1], n2[2], norm2);
+ double cosa;prosca (norm1, norm2, &cosa);
return cosa > min_cos;
}
@@ -188,67 +188,67 @@ bool evalSwapForOptimize(BDS_Edge *e, GFace *gf, BDS_Mesh &m)
{
if(e->numfaces() != 2) return false;
- BDS_Point *p11,*p12,*p13;
- BDS_Point *p21,*p22,*p23;
- BDS_Point *p31,*p32,*p33;
- BDS_Point *p41,*p42,*p43;
- swap_config(e,&p11,&p12,&p13,&p21,&p22,&p23,&p31,&p32,&p33,&p41,&p42,&p43);
+ BDS_Point *p11, *p12, *p13;
+ BDS_Point *p21, *p22, *p23;
+ BDS_Point *p31, *p32, *p33;
+ BDS_Point *p41, *p42, *p43;
+ swap_config(e, &p11, &p12, &p13, &p21, &p22, &p23, &p31, &p32, &p33, &p41, &p42, &p43);
// First, evaluate what we gain in element quality if the
// swap is performed
- double qa1 = qmTriangle(p11, p12, p13,QMTRI_RHO);
- double qa2 = qmTriangle(p21, p22, p23,QMTRI_RHO);
- double qb1 = qmTriangle(p31, p32, p33,QMTRI_RHO);
- double qb2 = qmTriangle(p41, p42, p43,QMTRI_RHO);
- double qa = std::min(qa1,qa2);
- double qb = std::min(qb1,qb2);
+ double qa1 = qmTriangle(p11, p12, p13, QMTRI_RHO);
+ double qa2 = qmTriangle(p21, p22, p23, QMTRI_RHO);
+ double qb1 = qmTriangle(p31, p32, p33, QMTRI_RHO);
+ double qb2 = qmTriangle(p41, p42, p43, QMTRI_RHO);
+ double qa = std::min(qa1, qa2);
+ double qb = std::min(qb1, qb2);
double qualIndicator = qb - qa;
bool qualShouldSwap = qb > 2*qa;
- bool qualCouldSwap = !(qb < qa*.5) && qb > .05;
+ bool qualCouldSwap = !(qb < qa * .5) && qb > .05;
// then evaluate if swap produces smoother surfaces
double norm11[3];
double norm12[3];
double norm21[3];
double norm22[3];
- normal_triangle(p11,p12,p13,norm11);
- normal_triangle(p21,p22,p23,norm12);
- normal_triangle(p31,p32,p33,norm21);
- normal_triangle(p41,p42,p43,norm22);
- double cosa;prosca (norm11,norm12,&cosa);
- double cosb;prosca (norm21,norm22,&cosb);
+ normal_triangle(p11, p12, p13, norm11);
+ normal_triangle(p21, p22, p23, norm12);
+ normal_triangle(p31, p32, p33, norm21);
+ normal_triangle(p41, p42, p43, norm22);
+ double cosa; prosca(norm11, norm12, &cosa);
+ double cosb; prosca(norm21, norm22, &cosb);
double smoothIndicator = cosb - cosa;
- bool smoothShouldSwap = (cosa < 0.1 && cosb > 0.3);
- bool smoothCouldSwap = !(cosb < cosa*.5);
+ bool smoothShouldSwap = (cosa < 0.1 && cosb > 0.3);
+ bool smoothCouldSwap = !(cosb < cosa * .5);
- double la = computeEdgeLinearLength ( p11 , p12 );
- double la_ = computeEdgeLinearLength ( p11 , p12 , gf , m.scalingU , m.scalingV);
- double lb = computeEdgeLinearLength ( p13 , p23 );
- double lb_ = computeEdgeLinearLength ( p13 , p23 , gf , m.scalingU , m.scalingV);
+ double la = computeEdgeLinearLength(p11, p12);
+ double la_ = computeEdgeLinearLength(p11, p12, gf, m.scalingU, m.scalingV);
+ double lb = computeEdgeLinearLength(p13, p23);
+ double lb_ = computeEdgeLinearLength(p13, p23, gf, m.scalingU, m.scalingV);
- double LA = (la_-la)/la_;
- double LB = (lb_-lb)/lb_;
+ double LA = (la_ -la) / la_;
+ double LB = (lb_ -lb) / lb_;
double distanceIndicator = LA - LB;
- bool distanceShouldSwap = (LB < .5*LA) && LA > 1.e-2;
- bool distanceCouldSwap = !(LB > 2*LA) || LB < 1.e-2;
+ bool distanceShouldSwap = (LB < .5 * LA) && LA > 1.e-2;
+ bool distanceCouldSwap = !(LB > 2 * LA) || LB < 1.e-2;
- if (20*qa < qb)return true;
-// if (LB > .025 && distanceIndicator > 0 && qb > .25)return true;
-// if (LB > .05 && distanceIndicator > 0 && qb > .15)return true;
-// if (LB > .1 && distanceIndicator > 0 && qb > .1)return true;
-// if (LB > .2 && distanceIndicator > 0 && qb > .05)return true;
-// if (LB > .3 && distanceIndicator > 0 && qb > .025)return true;
-
- // if swap enhances both criterion, the do it !
- if (distanceIndicator > 0 && qualIndicator > 0)return true;
- if (distanceShouldSwap && qualCouldSwap)return true;
- if (distanceCouldSwap && qualShouldSwap)return true;
-// if (smoothIndicator > 0 && qualIndicator > 0)return true;
-// if (smoothShouldSwap && qualCouldSwap)return true;
-// if (smoothCouldSwap && qualShouldSwap)return true;
- // if (distanceShouldSwap && qualCouldSwap)return true;
- // if (distanceCouldSwap && qualShouldSwap)return true;
+ if (20 * qa < qb) return true;
+ // if (LB > .025 && distanceIndicator > 0 && qb > .25) return true;
+ // if (LB > .05 && distanceIndicator > 0 && qb > .15) return true;
+ // if (LB > .1 && distanceIndicator > 0 && qb > .1) return true;
+ // if (LB > .2 && distanceIndicator > 0 && qb > .05) return true;
+ // if (LB > .3 && distanceIndicator > 0 && qb > .025) return true;
+
+ // if swap enhances both criterion, the do it!
+ if (distanceIndicator > 0 && qualIndicator > 0) return true;
+ if (distanceShouldSwap && qualCouldSwap) return true;
+ if (distanceCouldSwap && qualShouldSwap) return true;
+ // if (smoothIndicator > 0 && qualIndicator > 0) return true;
+ // if (smoothShouldSwap && qualCouldSwap) return true;
+ // if (smoothCouldSwap && qualShouldSwap) return true;
+ // if (distanceShouldSwap && qualCouldSwap) return true;
+ // if (distanceCouldSwap && qualShouldSwap) return true;
if (cosa < 0 && cosb > 0 && qb > 0.0)
return true;
return false;
@@ -256,28 +256,26 @@ bool evalSwapForOptimize(BDS_Edge *e, GFace *gf, BDS_Mesh &m)
bool edgeSwapTestDelaunay(BDS_Edge *e,GFace *gf)
{
-
BDS_Point *op[2];
- if(!e->p1->config_modified && ! e->p2->config_modified) return false;
+ if(!e->p1->config_modified && !e->p2->config_modified) return false;
if(e->numfaces() != 2) return false;
- e->oppositeof (op);
+ e->oppositeof(op);
- double p1x[3] = {e->p1->X,e->p1->Y,e->p1->Z};
- double p2x[3] = {e->p2->X,e->p2->Y,e->p2->Z};
- double op1x[3] = {op[0]->X,op[0]->Y,op[0]->Z};
- double op2x[3] = {op[1]->X,op[1]->Y,op[1]->Z};
+ double p1x[3] = {e->p1->X, e->p1->Y, e->p1->Z};
+ double p2x[3] = {e->p2->X, e->p2->Y, e->p2->Z};
+ double op1x[3] = {op[0]->X, op[0]->Y, op[0]->Z};
+ double op2x[3] = {op[1]->X, op[1]->Y, op[1]->Z};
double fourth[3];
- fourthPoint(p1x,p2x,op1x,fourth);
- double result = gmsh::insphere(p1x, p2x, op1x, fourth, op2x) * gmsh::orient3d(p1x, p2x, op1x, fourth);
- //double result = gmsh::incircle(p1u, p2u, op1u, op2u) * gmsh::orient2d(p1u, p2u, op1u);
- // printf("result = a%12.5E\n",result);
+ fourthPoint(p1x, p2x, op1x, fourth);
+ double result = gmsh::insphere(p1x, p2x, op1x, fourth, op2x) *
+ gmsh::orient3d(p1x, p2x, op1x, fourth);
return result > 0.;
}
-bool edgeSwapTestDelaunayAniso(BDS_Edge *e,GFace *gf,std::set<swapquad> & configs)
+bool edgeSwapTestDelaunayAniso(BDS_Edge *e, GFace *gf, std::set<swapquad> &configs)
{
BDS_Point *op[2];
@@ -287,43 +285,41 @@ bool edgeSwapTestDelaunayAniso(BDS_Edge *e,GFace *gf,std::set<swapquad> & config
e->oppositeof (op);
- swapquad sq (e->p1->iD,e->p2->iD,op[0]->iD,op[1]->iD);
- if (configs.find(sq) != configs.end())return false;
+ swapquad sq (e->p1->iD, e->p2->iD, op[0]->iD, op[1]->iD);
+ if (configs.find(sq) != configs.end()) return false;
configs.insert(sq);
- double edgeCenter[2] ={0.5*(e->p1->u+e->p2->u),
- 0.5*(e->p1->v+e->p2->v)};
+ double edgeCenter[2] ={0.5 * (e->p1->u + e->p2->u),
+ 0.5 * (e->p1->v + e->p2->v)};
- double p1[2] ={e->p1->u,e->p1->v};
- double p2[2] ={e->p2->u,e->p2->v};
- double p3[2] ={op[0]->u,op[0]->v};
- double p4[2] ={op[1]->u,op[1]->v};
+ double p1[2] ={e->p1->u, e->p1->v};
+ double p2[2] ={e->p2->u, e->p2->v};
+ double p3[2] ={op[0]->u, op[0]->v};
+ double p4[2] ={op[1]->u, op[1]->v};
double metric[3];
- buildMetric ( gf , edgeCenter , metric);
- // printf("%22.15E %22.15E %22.15E\n",metric[0],metric[1],metric[2]);
- if (!inCircumCircleAniso (gf,p1,p2,p3,p4,metric)){
+ buildMetric(gf, edgeCenter, metric);
+ if (!inCircumCircleAniso(gf, p1, p2, p3, p4, metric)){
return false;
}
return true;
}
-
bool evalSwap(BDS_Edge *e, double &qa, double &qb)
{
BDS_Point *op[2];
if(e->numfaces() != 2) return false;
e->oppositeof (op);
- double qa1 = qmTriangle(e->p1, e->p2, op[0],QMTRI_RHO);
- double qa2 = qmTriangle(e->p1, e->p2, op[1],QMTRI_RHO);
- double qb1 = qmTriangle(e->p1, op[0], op[1],QMTRI_RHO);
- double qb2 = qmTriangle(e->p2, op[0], op[1],QMTRI_RHO);
- qa = std::min(qa1,qa2);
- qb = std::min(qb1,qb2);
+ double qa1 = qmTriangle(e->p1, e->p2, op[0], QMTRI_RHO);
+ double qa2 = qmTriangle(e->p1, e->p2, op[1], QMTRI_RHO);
+ double qb1 = qmTriangle(e->p1, op[0], op[1], QMTRI_RHO);
+ double qb2 = qmTriangle(e->p2, op[0], op[1], QMTRI_RHO);
+ qa = std::min(qa1, qa2);
+ qb = std::min(qb1, qb2);
return true;
}
-int edgeSwapTestQuality(BDS_Edge *e, double fact = 1.1, bool force = false)
+int edgeSwapTestQuality(BDS_Edge *e, double fact=1.1, bool force=false)
{
BDS_Point *op[2];
@@ -335,30 +331,27 @@ int edgeSwapTestQuality(BDS_Edge *e, double fact = 1.1, bool force = false)
e->oppositeof (op);
if (!force)
- if (! edgeSwapTestAngle(e, cos(CTX.mesh.allow_swap_edge_angle*M_PI/180.)) ) return -1;
-
- double qa1 = qmTriangle(e->p1, e->p2, op[0],QMTRI_RHO);
- double qa2 = qmTriangle(e->p1, e->p2, op[1],QMTRI_RHO);
- double qb1 = qmTriangle(e->p1, op[0], op[1],QMTRI_RHO);
- double qb2 = qmTriangle(e->p2, op[0], op[1],QMTRI_RHO);
- double qa = std::min(qa1,qa2);
- double qb = std::min(qb1,qb2);
- if(qb > fact*qa) return 1;
- if(qb < qa/fact) return -1;
+ if (!edgeSwapTestAngle(e, cos(CTX.mesh.allow_swap_edge_angle * M_PI / 180.)))
+ return -1;
+
+ double qa1 = qmTriangle(e->p1, e->p2, op[0], QMTRI_RHO);
+ double qa2 = qmTriangle(e->p1, e->p2, op[1], QMTRI_RHO);
+ double qb1 = qmTriangle(e->p1, op[0], op[1], QMTRI_RHO);
+ double qb2 = qmTriangle(e->p2, op[0], op[1], QMTRI_RHO);
+ double qa = std::min(qa1, qa2);
+ double qb = std::min(qb1, qb2);
+ if(qb > fact * qa) return 1;
+ if(qb < qa / fact) return -1;
return 0;
}
-
-
-
-
-void swapEdgePass ( GFace *gf, BDS_Mesh &m, int &nb_swap )
+void swapEdgePass(GFace *gf, BDS_Mesh &m, int &nb_swap)
{
int NN1 = m.edges.size();
int NN2 = 0;
std::list<BDS_Edge*>::iterator it = m.edges.begin();
while (1){
- if (NN2++ >= NN1)break;
+ if (NN2++ >= NN1) break;
// result = -1 => forbid swap because too badly shaped elements
// result = 0 => whatever
// result = 1 => oblige to swap because the quality is greatly improved
@@ -366,18 +359,18 @@ void swapEdgePass ( GFace *gf, BDS_Mesh &m, int &nb_swap )
const double qual = CTX.mesh.algo2d == ALGO_2D_MESHADAPT ? 1 : 5;
int result = edgeSwapTestQuality(*it,qual);
if (CTX.mesh.algo2d == ALGO_2D_MESHADAPT )
- { if (m.swap_edge ( *it , BDS_SwapEdgeTestQuality(true)))nb_swap++; }
+ { if (m.swap_edge(*it, BDS_SwapEdgeTestQuality(true)))nb_swap++; }
else if ( result >= 0 && edgeSwapTestDelaunay(*it,gf))
- { if (m.swap_edge ( *it , BDS_SwapEdgeTestQuality(false))) nb_swap++; }
+ { if (m.swap_edge(*it, BDS_SwapEdgeTestQuality(false))) nb_swap++; }
}
++it;
}
}
-void gmshDelaunayizeBDS ( GFace *gf, BDS_Mesh &m, int &nb_swap )
+void gmshDelaunayizeBDS(GFace *gf, BDS_Mesh &m, int &nb_swap)
{
nb_swap = 0;
- std::set<swapquad> configs;
+ std::set<swapquad> configs;
while(1){
int NN1 = m.edges.size();
int NN2 = 0;
@@ -386,8 +379,8 @@ void gmshDelaunayizeBDS ( GFace *gf, BDS_Mesh &m, int &nb_swap )
while (1){
if (NN2++ >= NN1)break;
if (!(*it)->deleted){
- if (edgeSwapTestDelaunayAniso(*it,gf,configs)){
- if (m.swap_edge ( *it , BDS_SwapEdgeTestQuality(false))){
+ if (edgeSwapTestDelaunayAniso(*it, gf, configs)){
+ if (m.swap_edge(*it , BDS_SwapEdgeTestQuality(false))){
NSW++;
}
}
@@ -395,34 +388,34 @@ void gmshDelaunayizeBDS ( GFace *gf, BDS_Mesh &m, int &nb_swap )
++it;
}
nb_swap += NSW;
- if (!NSW)return;
+ if (!NSW) return;
}
}
-void splitEdgePassUnsorted ( GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
+void splitEdgePassUnsorted(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
{
int NN1 = m.edges.size();
int NN2 = 0;
std::list<BDS_Edge*>::iterator it = m.edges.begin();
while (1){
- if (NN2++ >= NN1)break;
+ if (NN2++ >= NN1) break;
if (!(*it)->deleted){
- double lone = NewGetLc ( *it,gf,m.scalingU,m.scalingV);
- if ((*it)->numfaces() == 2 && (lone > MAXE_)){
-
+ double lone = NewGetLc(*it, gf, m.scalingU, m.scalingV);
+ if ((*it)->numfaces() == 2 && (lone > MAXE_)){
//const double coord = 0.5;
- const double coord = computeEdgeMiddleCoord((*it)->p1,(*it)->p2,gf,m.scalingU,m.scalingV);
- BDS_Point *mid ;
+ const double coord = computeEdgeMiddleCoord((*it)->p1, (*it)->p2, gf,
+ m.scalingU, m.scalingV);
+ BDS_Point *mid;
mid = m.add_point(++m.MAXPOINTNUMBER,
coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u,
- coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v,gf);
-
- mid->lcBGM() = BGM_MeshSize(gf,
- (coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u)*m.scalingU,
- (coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v)*m.scalingV,
- mid->X,mid->Y,mid->Z);
- mid->lc() = 0.5 * ( (*it)->p1->lc() + (*it)->p2->lc() );
- if(!m.split_edge ( *it, mid )) m.del_point(mid);
+ coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v, gf);
+ mid->lcBGM() = BGM_MeshSize
+ (gf,
+ (coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u)*m.scalingU,
+ (coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v)*m.scalingV,
+ mid->X,mid->Y,mid->Z);
+ mid->lc() = 0.5 * ((*it)->p1->lc() + (*it)->p2->lc());
+ if(!m.split_edge(*it, mid)) m.del_point(mid);
else nb_split++;
}
}
@@ -430,45 +423,43 @@ void splitEdgePassUnsorted ( GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split
}
}
-void splitEdgePass ( GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
+void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
{
std::list<BDS_Edge*>::iterator it = m.edges.begin();
std::vector<std::pair<double, BDS_Edge*> > edges;
while (it != m.edges.end()){
if(!(*it)->deleted && (*it)->numfaces() == 2){
- double lone = NewGetLc(*it, gf,m.scalingU,m.scalingV);
+ double lone = NewGetLc(*it, gf, m.scalingU, m.scalingV);
if(lone > MAXE_){
- edges.push_back (std::make_pair (-lone, *it) );
+ edges.push_back(std::make_pair(-lone, *it));
}
}
++it;
}
- std::sort(edges.begin(),edges.end());
+ std::sort(edges.begin(), edges.end());
- for (int i=0;i<edges.size();++i){
+ for (int i = 0; i < edges.size(); ++i){
BDS_Edge *e = edges[i].second;
if (!e->deleted){
const double coord = 0.5;
- //const double coord = computeEdgeMiddleCoord(e->p1,e->p2,gf,m.scalingU,m.scalingV);
BDS_Point *mid ;
mid = m.add_point(++m.MAXPOINTNUMBER,
coord * e->p1->u + (1 - coord) * e->p2->u,
coord * e->p1->v + (1 - coord) * e->p2->v,gf);
-
- mid->lcBGM() = BGM_MeshSize(gf,
- (coord * e->p1->u + (1 - coord) * e->p2->u)*m.scalingU,
- (coord * e->p1->v + (1 - coord) * e->p2->v)*m.scalingV,
- mid->X,mid->Y,mid->Z);
- mid->lc() = 0.5 * ( e->p1->lc() + e->p2->lc() );
- if(!m.split_edge ( e, mid )) m.del_point(mid);
+ mid->lcBGM() = BGM_MeshSize
+ (gf,
+ (coord * e->p1->u + (1 - coord) * e->p2->u)*m.scalingU,
+ (coord * e->p1->v + (1 - coord) * e->p2->v)*m.scalingV,
+ mid->X,mid->Y,mid->Z);
+ mid->lc() = 0.5 * (e->p1->lc() + e->p2->lc());
+ if(!m.split_edge(e, mid)) m.del_point(mid);
else nb_split++;
}
}
}
-
void collapseEdgePass(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP, int &nb_collaps)
{
std::list<BDS_Edge*>::iterator it = m.edges.begin();
@@ -478,31 +469,30 @@ void collapseEdgePass(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP, int &nb_c
if(!(*it)->deleted && (*it)->numfaces() == 2){
double lone = NewGetLc(*it, gf,m.scalingU,m.scalingV);
if(lone < MINE_){
- edges.push_back (std::make_pair (lone, *it) );
+ edges.push_back (std::make_pair(lone, *it));
}
}
++it;
}
- std::sort(edges.begin(),edges.end());
+ std::sort(edges.begin(), edges.end());
- for (int i=0;i<edges.size();i++){
+ for (int i = 0; i < edges.size(); i++){
BDS_Edge *e = edges[i].second;
- // printf("%12.5E\n",edges[i].first);
if(!e->deleted){
bool res = false;
if(e->p1->iD > MAXNP)
- res = m.collapse_edge_parametric(e,e->p1);
+ res = m.collapse_edge_parametric(e, e->p1);
else if(e->p2->iD > MAXNP)
- res = m.collapse_edge_parametric(e,e->p2);
+ res = m.collapse_edge_parametric(e, e->p2);
if(res)
nb_collaps++;
}
}
}
-
-void collapseEdgePassUnSorted(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP, int &nb_collaps)
+void collapseEdgePassUnSorted(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP,
+ int &nb_collaps)
{
int NN1 = m.edges.size();
int NN2 = 0;
@@ -512,7 +502,7 @@ void collapseEdgePassUnSorted(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP, i
if(NN2++ >= NN1) break;
if(!(*it)->deleted){
- double lone = NewGetLc(*it, gf,m.scalingU,m.scalingV);
+ double lone = NewGetLc(*it, gf, m.scalingU, m.scalingV);
if(!(*it)->deleted && (*it)->numfaces() == 2 && lone < MINE_){
bool res = false;
if((*it)->p1->iD > MAXNP)
@@ -527,10 +517,8 @@ void collapseEdgePassUnSorted(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP, i
}
}
-
void smoothVertexPass(GFace *gf, BDS_Mesh &m, int &nb_smooth, bool q)
{
- // return;
std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin();
while(itp != m.points.end()){
if(m.smooth_point_centroid(*itp, gf,q))
@@ -539,13 +527,11 @@ void smoothVertexPass(GFace *gf, BDS_Mesh &m, int &nb_smooth, bool q)
}
}
-void gmshRefineMeshBDS (GFace *gf,
- BDS_Mesh &m,
- const int NIT,
- const bool computeNodalSizeField)
+void gmshRefineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
+ const bool computeNodalSizeField)
{
int IT = 0;
-
+
int MAXNP = m.MAXPOINTNUMBER;
// IF ASKED , compute nodal size field using 1D Mesh
@@ -554,12 +540,12 @@ void gmshRefineMeshBDS (GFace *gf,
while (itp != m.points.end()){
std::list<BDS_Edge*>::iterator it = (*itp)->edges.begin();
std::list<BDS_Edge*>::iterator ite = (*itp)->edges.end();
- double L=0;
+ double L = 0;
int ne = 0;
- while(it!=ite){
+ while(it != ite){
double l = (*it)->length();
if ((*it)->g && (*it)->g->classif_degree == 1){
- L=ne?std::max(L,l):l;
+ L = ne ? std::max(L, l) : l;
ne++;
}
++it;
@@ -572,99 +558,88 @@ void gmshRefineMeshBDS (GFace *gf,
}
}
- double t_spl=0, t_sw=0,t_col=0,t_sm=0;
+ double t_spl = 0, t_sw = 0,t_col = 0,t_sm = 0;
- const double MINE_ = 0.67, MAXE_=1.4;
+ const double MINE_ = 0.67, MAXE_ = 1.4;
double mesh_quality = 0;
- while (1)
- {
- // we count the number of local mesh modifs.
-
- int nb_split=0;
- int nb_smooth=0;
- int nb_collaps=0;
- int nb_swap=0;
-
- // split long edges
-
- double minL=1.E22,maxL=0;
-
- int NN1 = m.edges.size();
- int NN2 = 0;
- std::list<BDS_Edge*>::iterator it = m.edges.begin();
- while (1){
- if (NN2++ >= NN1)break;
- if (!(*it)->deleted)
- {
- (*it)->p1->config_modified = false;
- (*it)->p2->config_modified = false;
- double lone = NewGetLc ( *it,gf,m.scalingU,m.scalingV);
- maxL = std::max(maxL,lone);
- minL = std::min(minL,lone);
- }
- ++it;
+ while (1){
+ // we count the number of local mesh modifs.
+ int nb_split =0;
+ int nb_smooth = 0 ;
+ int nb_collaps =0;
+ int nb_swap = 0;
+
+ // split long edges
+ double minL = 1.E22, maxL = 0;
+ int NN1 = m.edges.size();
+ int NN2 = 0;
+ std::list<BDS_Edge*>::iterator it = m.edges.begin();
+ while (1){
+ if (NN2++ >= NN1)break;
+ if (!(*it)->deleted){
+ (*it)->p1->config_modified = false;
+ (*it)->p2->config_modified = false;
+ double lone = NewGetLc(*it, gf, m.scalingU, m.scalingV);
+ maxL = std::max(maxL, lone);
+ minL = std::min(minL, lone);
}
-
- if ((minL > MINE_ && maxL < MAXE_) || IT > (abs(NIT)))break;
- double maxE = MAXE_;
- double minE = MINE_;
- double t1 = Cpu();
- splitEdgePass ( gf, m, maxE, nb_split);
- double t2 = Cpu();
- swapEdgePass ( gf, m, nb_swap);
- swapEdgePass ( gf, m, nb_swap);
- swapEdgePass ( gf, m, nb_swap);
- double t3 = Cpu();
- collapseEdgePass ( gf, m, minE , MAXNP, nb_collaps);
- double t4 = Cpu();
- // swapEdgePass ( gf, m, nb_swap);
- double t5 = Cpu();
- smoothVertexPass ( gf, m, nb_smooth,false);
- double t6 = Cpu();
- // swapEdgePass ( gf, m, nb_swap);
- double t7 = Cpu();
- // clean up the mesh
-
- t_spl += t2 - t1;
- t_sw += t3 - t2;
- t_sw += t5 - t4;
- t_sw += t7 - t6;
- t_col += t4 - t3;
- t_sm += t6 - t5;
-
- m.cleanup();
- IT++;
- Msg(DEBUG1," iter %3d minL %8.3f/%8.3f maxL %8.3f/%8.3f : %6d splits, %6d swaps, %6d collapses, %6d moves",IT,minL,minE,maxL,maxE,nb_split,nb_swap,nb_collaps,nb_smooth);
-
- if (nb_split==0 && nb_collaps == 0)break;
- }
-
+ ++it;
+ }
+ if ((minL > MINE_ && maxL < MAXE_) || IT > (abs(NIT))) break;
+ double maxE = MAXE_;
+ double minE = MINE_;
+ double t1 = Cpu();
+ splitEdgePass(gf, m, maxE, nb_split);
+ double t2 = Cpu();
+ swapEdgePass(gf, m, nb_swap);
+ swapEdgePass(gf, m, nb_swap);
+ swapEdgePass(gf, m, nb_swap);
+ double t3 = Cpu();
+ collapseEdgePass(gf, m, minE, MAXNP, nb_collaps);
+ double t4 = Cpu();
+ double t5 = Cpu();
+ smoothVertexPass(gf, m, nb_smooth, false);
+ double t6 = Cpu();
+ // swapEdgePass ( gf, m, nb_swap);
+ double t7 = Cpu();
+ // clean up the mesh
+ t_spl += t2 - t1;
+ t_sw += t3 - t2;
+ t_sw += t5 - t4;
+ t_sw += t7 - t6;
+ t_col += t4 - t3;
+ t_sm += t6 - t5;
+ m.cleanup();
+ IT++;
+ Msg(DEBUG1," iter %3d minL %8.3f/%8.3f maxL %8.3f/%8.3f : "
+ "%6d splits, %6d swaps, %6d collapses, %6d moves",
+ IT, minL, minE, maxL, maxE, nb_split, nb_swap, nb_collaps, nb_smooth);
+ if (nb_split == 0 && nb_collaps == 0) break;
+ }
+
double t_total = t_spl + t_sw + t_col + t_sm;
if(!t_total) t_total = 1.e-6;
Msg(DEBUG1," ---------------------------------------");
Msg(DEBUG1," CPU Report ");
Msg(DEBUG1," ---------------------------------------");
- Msg(DEBUG1," CPU SWAP %12.5E sec (%3.0f %%)",t_sw,100 * t_sw/t_total);
- Msg(DEBUG1," CPU SPLIT %12.5E sec (%3.0f %%) ",t_spl,100 * t_spl/t_total);
- Msg(DEBUG1," CPU COLLAPS %12.5E sec (%3.0f %%) ",t_col,100 * t_col/t_total);
- Msg(DEBUG1," CPU SMOOTH %12.5E sec (%3.0f %%) ",t_sm,100 * t_sm/t_total);
+ Msg(DEBUG1," CPU SWAP %12.5E sec (%3.0f %%)", t_sw,100 * t_sw / t_total);
+ Msg(DEBUG1," CPU SPLIT %12.5E sec (%3.0f %%) ", t_spl,100 * t_spl / t_total);
+ Msg(DEBUG1," CPU COLLAPS %12.5E sec (%3.0f %%) ", t_col,100 * t_col / t_total);
+ Msg(DEBUG1," CPU SMOOTH %12.5E sec (%3.0f %%) ", t_sm,100 * t_sm / t_total);
Msg(DEBUG1," ---------------------------------------");
Msg(DEBUG1," CPU TOTAL %12.5E sec ",t_total);
Msg(DEBUG1," ---------------------------------------");
}
-void allowAppearanceofEdge (BDS_Point *p1, BDS_Point *p2){
+void allowAppearanceofEdge (BDS_Point *p1, BDS_Point *p2)
+{
}
-
-void invalidEdgesPeriodic(BDS_Mesh &m,
- std::map<BDS_Point*,MVertex*> *recover_map,
+void invalidEdgesPeriodic(BDS_Mesh &m, std::map<BDS_Point*,MVertex*> *recover_map,
std::set<BDS_Edge*> &toSplit)
{
-
// first look for degenerated vertices
-
std::list<BDS_Edge*>::iterator it = m.edges.begin();
std::set<MVertex *> degenerated;
while (it != m.edges.end()){
@@ -681,8 +656,6 @@ void invalidEdgesPeriodic(BDS_Mesh &m,
++it;
}
- // printf("%d degenerated\n",degenerated.size());
-
toSplit.clear();
it = m.edges.begin();
std::map< std::pair < MVertex*, BDS_Point*> , BDS_Edge *> touchPeriodic;
@@ -696,22 +669,25 @@ void invalidEdgesPeriodic(BDS_Mesh &m,
itp1->second == itp2->second) toSplit.insert(e);
else if (itp1 != recover_map->end() && itp2 == recover_map->end()){
std::pair<MVertex*,BDS_Point*> a ( itp1->second, e->p2 );
- std::map< std::pair < MVertex*, BDS_Point*> , BDS_Edge *> :: iterator itf =
- touchPeriodic.find (a);
+ std::map<std::pair<MVertex*, BDS_Point*>, BDS_Edge*>::iterator itf =
+ touchPeriodic.find(a);
if (itf == touchPeriodic.end()) touchPeriodic[a] = e;
- else if (degenerated.find(itp1->second)==degenerated.end()){toSplit.insert(e);toSplit.insert(itf->second);}
+ else if (degenerated.find(itp1->second) == degenerated.end()){
+ toSplit.insert(e); toSplit.insert(itf->second);
+ }
}
else if (itp1 == recover_map->end() && itp2 != recover_map->end()){
- std::pair<MVertex*,BDS_Point*> a ( itp2->second, e->p1 );
- std::map< std::pair < MVertex*, BDS_Point*> , BDS_Edge *> :: iterator itf =
+ std::pair<MVertex*,BDS_Point*> a(itp2->second, e->p1);
+ std::map<std::pair<MVertex*, BDS_Point*>, BDS_Edge*>::iterator itf =
touchPeriodic.find (a);
if (itf == touchPeriodic.end()) touchPeriodic[a] = e;
- else if (degenerated.find(itp2->second)==degenerated.end()){toSplit.insert(e);toSplit.insert(itf->second);}
+ else if (degenerated.find(itp2->second) == degenerated.end()){
+ toSplit.insert(e); toSplit.insert(itf->second);
+ }
}
}
++it;
}
- // printf("%d edges to splitounette\n",toSplit.size());
}
// consider p1 and p2, 2 vertices that are different in the parametric
@@ -721,128 +697,130 @@ void invalidEdgesPeriodic(BDS_Mesh &m,
// if not do not allow p1 v and p2 v, split them
// if p1 p2 exists and it is internal, split it
-int gmshSolveInvalidPeriodic(GFace *gf,
- BDS_Mesh &m,
- std::map<BDS_Point*,MVertex*> *recover_map){
- // return 0;
+int gmshSolveInvalidPeriodic(GFace *gf, BDS_Mesh &m,
+ std::map<BDS_Point*,MVertex*> *recover_map)
+{
std::set<BDS_Edge*> toSplit;
- invalidEdgesPeriodic(m,recover_map,toSplit);
+ invalidEdgesPeriodic(m, recover_map, toSplit);
std::set<BDS_Edge*>::iterator ite = toSplit.begin();
- // printf("%d edges to split\n",toSplit.size());
- for (;ite !=toSplit.end();++ite){
+
+ for (;ite !=toSplit.end(); ++ite){
BDS_Edge *e = *ite;
if (!e->deleted && e->numfaces() == 2){
const double coord = 0.5;
BDS_Point *mid ;
- // printf("%d %d\n",e->p1->iD,e->p2->iD);
- mid = m.add_point(++m.MAXPOINTNUMBER,
- coord * e->p1->u + (1 - coord) * e->p2->u,
- coord * e->p1->v + (1 - coord) * e->p2->v,gf);
+ mid = m.add_point(++m.MAXPOINTNUMBER,
+ coord * e->p1->u + (1 - coord) * e->p2->u,
+ coord * e->p1->v + (1 - coord) * e->p2->v, gf);
mid->lcBGM() = BGM_MeshSize(gf,
- (coord * e->p1->u + (1 - coord) * e->p2->u)*m.scalingU,
- (coord * e->p1->v + (1 - coord) * e->p2->v)*m.scalingV,
- mid->X,mid->Y,mid->Z);
- mid->lc() = 0.5 * ( e->p1->lc() + e->p2->lc() );
- if(!m.split_edge ( e, mid )) m.del_point(mid);
+ (coord * e->p1->u + (1 - coord) * e->p2->u) * m.scalingU,
+ (coord * e->p1->v + (1 - coord) * e->p2->v) * m.scalingV,
+ mid->X, mid->Y, mid->Z);
+ mid->lc() = 0.5 * (e->p1->lc() + e->p2->lc());
+ if(!m.split_edge(e, mid)) m.del_point(mid);
}
}
int nb_smooth;
- if (toSplit.size())smoothVertexPass ( gf,m,nb_smooth,true);
- m.cleanup();
+ if(toSplit.size()) smoothVertexPass(gf, m, nb_smooth, true);
+ m.cleanup();
return toSplit.size();
}
-void gmshOptimizeMeshBDS(GFace *gf,
- BDS_Mesh &m,
- const int NIT,
- std::map<BDS_Point*,MVertex*> *recover_map = 0)
+void gmshOptimizeMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
+ std::map<BDS_Point*,MVertex*> *recover_map=0)
{
int nb_swap;
- gmshDelaunayizeBDS ( gf, m, nb_swap );
+ gmshDelaunayizeBDS(gf, m, nb_swap);
- for (int ITER = 0;ITER < 3;ITER++){
+ for (int ITER = 0; ITER < 3; ITER++){
double LIMIT = .1;
- for (int KK=0;KK<4;KK++){
+ for (int KK = 0; KK < 4; KK++){
// swap edges that provide a better configuration
int NN1 = m.edges.size();
int NN2 = 0;
std::list<BDS_Edge*>::iterator it = m.edges.begin();
while (1){
if (NN2++ >= NN1)break;
- if (evalSwapForOptimize(*it,gf,m))
- m.swap_edge ( *it , BDS_SwapEdgeTestQuality(false));
+ if (evalSwapForOptimize(*it, gf, m))
+ m.swap_edge(*it, BDS_SwapEdgeTestQuality(false));
++it;
}
m.cleanup();
int nb_smooth;
- smoothVertexPass ( gf,m,nb_smooth,true);
+ smoothVertexPass(gf, m, nb_smooth, true);
}
}
-
+
int nbSplit = 0;
if (recover_map){
- while(gmshSolveInvalidPeriodic(gf,m,recover_map)){}
+ while(gmshSolveInvalidPeriodic(gf, m, recover_map)){
+ }
}
}
// DELAUNAY BDS
-void delaunayPointInsertionBDS ( GFace *gf, BDS_Mesh &m, BDS_Point *v, BDS_Face *f){
- const double p[2] = {v->u,v->v};
+void delaunayPointInsertionBDS(GFace *gf, BDS_Mesh &m, BDS_Point *v, BDS_Face *f)
+{
+ const double p[2] = {v->u, v->v};
BDS_Edge *e1 = f->e1;
BDS_Edge *e2 = f->e2;
BDS_Edge *e3 = f->e3;
// face is splitted,
- m.split_face ( f , v );
+ m.split_face(f, v);
int nb_swap = 0;
- gmshDelaunayizeBDS ( gf, m, nb_swap );
+ gmshDelaunayizeBDS(gf, m, nb_swap);
}
// build the BDS from a list of GFace
// This is a TRUE copy
-BDS_Mesh * gmsh2BDS ( std::list<GFace*> & l){
+BDS_Mesh *gmsh2BDS(std::list<GFace*> &l)
+{
BDS_Mesh *m = new BDS_Mesh;
- for (std::list<GFace*>::iterator it = l.begin();it!=l.end();++it){
+ for (std::list<GFace*>::iterator it = l.begin(); it != l.end(); ++it){
GFace *gf = *it;
- m->add_geom (gf->tag(), 2);
+ m->add_geom(gf->tag(), 2);
BDS_GeomEntity *g2 = m->get_geom(gf->tag(), 2);
- for (int i=0;i<gf->triangles.size();i++){
+ for (int i = 0; i < gf->triangles.size(); i++){
MTriangle *e = gf->triangles[i];
BDS_Point *p[3];
- for (int j=0;j<3;j++){
+ for (int j = 0; j < 3; j++){
p[j] = m->find_point(e->getVertex(j)->getNum());
if (!p[j]) {
- p[j] = m->add_point(e->getVertex(j)->getNum(),e->getVertex(j)->x(),e->getVertex(j)->y(),e->getVertex(j)->z());
- double u0,v0;
- parametricCoordinates ( e->getVertex(j), gf, u0, v0);
+ p[j] = m->add_point(e->getVertex(j)->getNum(), e->getVertex(j)->x(),
+ e->getVertex(j)->y(), e->getVertex(j)->z());
+ double u0, v0;
+ parametricCoordinates(e->getVertex(j), gf, u0, v0);
p[j]->u = u0;
p[j]->v = v0;
- m->add_geom (e->getVertex(j)->onWhat()->tag(), e->getVertex(j)->onWhat()->dim());
- BDS_GeomEntity *g = m->get_geom(e->getVertex(j)->onWhat()->tag(), e->getVertex(j)->onWhat()->dim());
+ m->add_geom(e->getVertex(j)->onWhat()->tag(),
+ e->getVertex(j)->onWhat()->dim());
+ BDS_GeomEntity *g = m->get_geom(e->getVertex(j)->onWhat()->tag(),
+ e->getVertex(j)->onWhat()->dim());
p[j]->g = g;
}
}
- BDS_Face *f = m->add_triangle ( p[0]->iD,p[1]->iD,p[2]->iD);
+ BDS_Face *f = m->add_triangle(p[0]->iD, p[1]->iD, p[2]->iD);
f->g = g2;
}
}
return m;
}
-void gmshCollapseSmallEdges (GModel &gm){
+void gmshCollapseSmallEdges(GModel &gm)
+{
return;
gm.renumberMeshVertices(true);
std::list<GFace*> faces;
for (GModel::fiter fit = gm.firstFace(); fit != gm.lastFace(); fit++){
faces.push_back(*fit);
}
- BDS_Mesh *pm = gmsh2BDS (faces);
- outputScalarField(pm->triangles,"all.pos",0);
+ BDS_Mesh *pm = gmsh2BDS(faces);
+ outputScalarField(pm->triangles, "all.pos", 0);
for (GModel::eiter eit = gm.firstEdge(); eit != gm.lastEdge(); eit++){
}
-
delete pm;
}
diff --git a/Mesh/meshGFaceBDS.h b/Mesh/meshGFaceBDS.h
index e5b5b470fa9fc91828bce0e8914a495d357a730d..ce181b1b691d6651654e3bac6d38573130518b2f 100644
--- a/Mesh/meshGFaceBDS.h
+++ b/Mesh/meshGFaceBDS.h
@@ -27,19 +27,15 @@ class BDS_Mesh;
class BDS_Point;
class MVertex;
-void computeMeshSizeFieldAccuracy (GFace *gf, BDS_Mesh &m, double &avg, double &max_e, double &min_e, int &nE, int &GS);
-void computeElementShapes (GFace *gf, BDS_Mesh &m, double &worst, double &avg, double &best, int &nT, int &nbGQ);
-void gmshRefineMeshBDS (GFace *gf,
- BDS_Mesh &m,
- const int NIT,
- const bool computeNodalSizeField);
-void gmshOptimizeMeshBDS(GFace *gf,
- BDS_Mesh &m,
- const int NIT,
- std::map<BDS_Point*,MVertex*> *recover_map = 0);
-
-void gmshDelaunayizeBDS ( GFace *gf, BDS_Mesh &m, int &nb_swap );
-
-void gmshCollapseSmallEdges (GModel &gm);
+void computeMeshSizeFieldAccuracy(GFace *gf, BDS_Mesh &m, double &avg,
+ double &max_e, double &min_e, int &nE, int &GS);
+void computeElementShapes(GFace *gf, BDS_Mesh &m, double &worst, double &avg,
+ double &best, int &nT, int &nbGQ);
+void gmshRefineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
+ const bool computeNodalSizeField);
+void gmshOptimizeMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
+ std::map<BDS_Point*,MVertex*> *recover_map=0);
+void gmshDelaunayizeBDS(GFace *gf, BDS_Mesh &m, int &nb_swap);
+void gmshCollapseSmallEdges(GModel &gm);
#endif
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index d19c750984d3921af7dc896507a164b77401d8d6..93d4516147614e9e78b8f7c814cde0061ba78d8d 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegionDelaunayInsertion.cpp,v 1.36 2008-02-21 09:49:22 geuzaine Exp $
+// $Id: meshGRegionDelaunayInsertion.cpp,v 1.37 2008-02-21 12:11:12 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -723,8 +723,8 @@ void gmshOptimizeMesh (GRegion *gr, const gmshQualityMeasure4Tet &qm)
void insertVerticesInRegion (GRegion *gr)
{
- printf("sizeof MTet4 = %d sizeof MTetrahedron %d sizeof(MVertex) %d\n",
- sizeof(MTet4),sizeof(MTetrahedron), sizeof(MVertex));
+ //printf("sizeof MTet4 = %d sizeof MTetrahedron %d sizeof(MVertex) %d\n",
+ // sizeof(MTet4),sizeof(MTetrahedron), sizeof(MVertex));
std::map<MVertex*,double> vSizesMap;
std::vector<double> vSizes;
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index 68aacd7806d50257b76438bc6a4879d9fffb8803..3901a702d74fe51537b80d96e3b58132947a7b02 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -1,4 +1,4 @@
-// $Id: qualityMeasures.cpp,v 1.8 2008-02-21 09:45:15 remacle Exp $
+// $Id: qualityMeasures.cpp,v 1.9 2008-02-21 12:11:12 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -71,9 +71,9 @@ double qmTriangle(const double &xa, const double &ya, const double &za,
const double lb = norme (b);
const double lc = norme (c);
- double pva [3]; prodve (b,c,pva); const double sina = norm3 (pva);
- double pvb [3]; prodve (c,a,pvb); const double sinb = norm3 (pvb);
- double pvc [3]; prodve (a,b,pvc); const double sinc = norm3 (pvc);
+ double pva [3]; prodve(b,c,pva); const double sina = norm3(pva);
+ double pvb [3]; prodve(c,a,pvb); const double sinb = norm3(pvb);
+ double pvc [3]; prodve(a,b,pvc); const double sinc = norm3(pvc);
if (sina == 0.0 && sinb == 0.0 && sinc == 0.0) quality = 0.0;
else quality = 2 * (2*sina*sinb*sinc/(sina + sinb + sinc) );
@@ -85,9 +85,9 @@ double qmTriangle(const double &xa, const double &ya, const double &za,
double a [3] = {xc-xa,yc-ya,zc-za};
double b [3] = {xb-xa,yb-ya,zb-za};
double c [3] ; prodve(a,b,c); norme(c);
- double A[3][3] = { a[0] , b[0] , c[0] ,
- a[1] , b[1] , c[1] ,
- a[2] , b[2] , c[2] };
+ double A[3][3] = {a[0] , b[0] , c[0] ,
+ a[1] , b[1] , c[1] ,
+ a[2] , b[2] , c[2]};
quality = -1;
}
break;
diff --git a/Numeric/FunctionSpace.cpp b/Numeric/FunctionSpace.cpp
index 9dc7527aa675b783612cfa0e101f4948a2c49341..b0dc408b15e8dd03ec26057f4b2f8935706ea2ce 100644
--- a/Numeric/FunctionSpace.cpp
+++ b/Numeric/FunctionSpace.cpp
@@ -1,83 +1,98 @@
+// $Id: FunctionSpace.cpp,v 1.2 2008-02-21 12:11:12 geuzaine Exp $
+//
+// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
#include <cmath>
#include "FunctionSpace.h"
#include "GmshDefines.h"
-Double_Matrix generatePascalTriangle(int order) {
-
- Double_Matrix monomials((order+1)*(order+2)/2,2);
+Double_Matrix generatePascalTriangle(int order)
+{
+ Double_Matrix monomials((order + 1) * (order + 2) / 2, 2);
int index = 0;
- for (int i=0;i<=order;i++) {
- for (int j=0;j<=i;j++) {
- monomials(index,0) = i - j;
- monomials(index,1) = j;
+ for (int i = 0; i <= order; i++) {
+ for (int j = 0; j <= i; j++) {
+ monomials(index, 0) = i - j;
+ monomials(index, 1) = j;
index++;
}
}
return monomials;
}
-// -----------------------------------------------------------------------------
-/*! generate the exterior hull of the Pascal triangle for Serendipity element */
-// -----------------------------------------------------------------------------
-
-Double_Matrix generatePascalSerendipityTriangle(int order) {
+// generate the exterior hull of the Pascal triangle for Serendipity element
- Double_Matrix monomials(3 * order,2);
+Double_Matrix generatePascalSerendipityTriangle(int order)
+{
+ Double_Matrix monomials(3 * order, 2);
- monomials(0,0) = 0;
- monomials(0,1) = 0;
+ monomials(0, 0) = 0;
+ monomials(0, 1) = 0;
int index = 1;
- for (int i=1;i<=order;i++) {
- if (i==order) {
- for (int j=0;j<=i;j++) {
- monomials(index,0) = i - j;
- monomials(index,1) = j;
+ for (int i = 1; i <= order; i++) {
+ if (i == order) {
+ for (int j = 0; j <= i; j++) {
+ monomials(index, 0) = i - j;
+ monomials(index, 1) = j;
index++;
}
}
else {
- monomials(index,0) = i;
- monomials(index,1) = 0;
+ monomials(index, 0) = i;
+ monomials(index, 1) = 0;
index++;
-
- monomials(index,0) = 0;
- monomials(index,1) = i;
+ monomials(index, 0) = 0;
+ monomials(index, 1) = i;
index++;
}
}
return monomials;
}
+// generate the monomials subspace of all monomials of order exactly == p
-// -----------------------------------------------------------------------------
-/*!\brief generate the monomials subspace of all monomials of order exactly == p */
-// -----------------------------------------------------------------------------
-
-Double_Matrix generateMonomialSubspace(int dim,int p) {
-
+Double_Matrix generateMonomialSubspace(int dim, int p)
+{
Double_Matrix monomials;
switch (dim) {
case 1:
- monomials = Double_Matrix(1,1);
- monomials(0,0) = p;
+ monomials = Double_Matrix(1, 1);
+ monomials(0, 0) = p;
break;
case 2:
- monomials = Double_Matrix(p+1,2);
- for (int k=0;k<=p;k++) {
- monomials(k,0) = p - k;
- monomials(k,1) = k;
+ monomials = Double_Matrix(p + 1, 2);
+ for (int k = 0; k <= p; k++) {
+ monomials(k, 0) = p - k;
+ monomials(k, 1) = k;
}
break;
case 3:
- monomials = Double_Matrix((p+1)*(p+2)/2,3);
+ monomials = Double_Matrix((p + 1) * (p + 2) / 2, 3);
int index = 0;
- for (int i=0;i<=p;i++) {
- for (int k=0;k<=p-i;k++) {
- monomials(index,0) = p - i - k;
- monomials(index,1) = k;
- monomials(index,2) = i;
+ for (int i = 0; i <= p; i++) {
+ for (int k = 0; k <= p - i; k++) {
+ monomials(index, 0) = p - i - k;
+ monomials(index, 1) = k;
+ monomials(index, 2) = i;
index++;
}
}
@@ -86,310 +101,332 @@ Double_Matrix generateMonomialSubspace(int dim,int p) {
return monomials;
}
+// generate external hull of the Pascal tetrahedron
-// -----------------------------------------------------------------------------
-/*!\brief generate external hull of the Pascal tetrahedron */
-// -----------------------------------------------------------------------------
-
-Double_Matrix generatePascalSerendipityTetrahedron(int order) {
-
- int nbMonomials = 4 + 6 * std::max(0,order-1) + 4 * std::max(0,(order-2)*(order-1)/2);
- Double_Matrix monomials(nbMonomials,3);
+Double_Matrix generatePascalSerendipityTetrahedron(int order)
+{
+ int nbMonomials = 4 + 6 * std::max(0, order - 1) +
+ 4 * std::max(0, (order - 2) * (order - 1) / 2);
+ Double_Matrix monomials(nbMonomials, 3);
// order 0
monomials.set_all(0);
-
- monomials(0,0) = 0;
- monomials(0,1) = 0;
- monomials(0,2) = 0;
+
+ monomials(0, 0) = 0;
+ monomials(0, 1) = 0;
+ monomials(0, 2) = 0;
int index = 1;
- for (int p=1;p<order;p++) {
- for (int i=0;i<3;i++) {
- int j = (i+1)%3;
- int k = (i+2)%3;
- for (int ii=0;ii<p;ii++) {
- monomials(index,i) = p - ii;
- monomials(index,j) = ii;
- monomials(index,k) = 0;
+ for (int p = 1; p < order; p++) {
+ for (int i = 0; i < 3; i++) {
+ int j = (i + 1) % 3;
+ int k = (i + 2) % 3;
+ for (int ii = 0; ii < p; ii++) {
+ monomials(index, i) = p - ii;
+ monomials(index, j) = ii;
+ monomials(index, k) = 0;
index++;
}
}
}
- Double_Matrix monomialsMaxOrder = generateMonomialSubspace(3,order);
+ Double_Matrix monomialsMaxOrder = generateMonomialSubspace(3, order);
int nbMaxOrder = monomialsMaxOrder.size1();
- Double_Matrix(monomials.touchSubmatrix(index,nbMaxOrder,0,3)).memcpy(monomialsMaxOrder);
+ Double_Matrix(monomials.touchSubmatrix(index, nbMaxOrder, 0, 3)).memcpy(monomialsMaxOrder);
return monomials;
}
-// -----------------------------------------------------------------------------
-/*!\brief generate Pascal tetrahedron */
-// -----------------------------------------------------------------------------
+// generate Pascal tetrahedron
+Double_Matrix generatePascalTetrahedron(int order)
+{
+ int nbMonomials = (order + 1) * (order + 2) * (order + 3) / 6;
-Double_Matrix generatePascalTetrahedron(int order) {
-
- int nbMonomials = (order+1)*(order+2)*(order+3) / 6;
-
- Double_Matrix monomials(nbMonomials,3);
+ Double_Matrix monomials(nbMonomials, 3);
int index = 0;
- for (int p=0;p<=order;p++) {
- Double_Matrix monOrder = generateMonomialSubspace(3,p);
+ for (int p = 0; p <= order; p++) {
+ Double_Matrix monOrder = generateMonomialSubspace(3, p);
int nb = monOrder.size1();
- Double_Matrix(monomials.touchSubmatrix(index,nb,0,3)).memcpy(monOrder);
+ Double_Matrix(monomials.touchSubmatrix(index, nb, 0, 3)).memcpy(monOrder);
index += nb;
}
return monomials;
}
+int nbdoftriangle(int order) { return (order + 1) * (order + 2) / 2; }
+int nbdoftriangleserendip(int order) { return 3 * order; }
-
-int nbdoftriangle(int order) {return (order+1)*(order+2)/2;}
-int nbdoftriangleserendip(int order) {return 3 * order;}
-
-void nodepositionface0(int order, double *u, double *v, double *w){
-
+void nodepositionface0(int order, double *u, double *v, double *w)
+{
int ndofT = nbdoftriangle(order);
- if (order ==0) {u[0]=0.;v[0]=0.; return;}
+ if (order == 0) { u[0] = 0.; v[0] = 0.; return; }
- u[0]= 0.; v[0]= 0.; w[0] = 0.;
- u[1]= order; v[1]= 0; w[1] = 0.;
+ u[0]= 0.; v[0]= 0.; w[0] = 0.;
+ u[1]= order; v[1]= 0; w[1] = 0.;
u[2]= 0.; v[2]= order; w[2] = 0.;
+
// edges
- for (int k=0; k< (order-1); k++){
- u[3+k] = k+1; v[3+k]=0.; w[3+k] = 0.;
- u[3+order-1+k] = order -1 -k ; v [3+order-1+k] = k+1; w[3+order-1+k] = 0.;
- u[3+2*(order-1)+k] = 0. ; v [3+2*(order-1)+k] = order -1 -k; w [3+2*(order-1)+k] = 0.;
+ for (int k = 0; k < (order - 1); k++){
+ u[3 + k] = k + 1;
+ v[3 + k] =0.;
+ w[3 + k] = 0.;
+
+ u[3 + order - 1 + k] = order - 1 - k ;
+ v[3 + order - 1 + k] = k + 1;
+ w[3 + order - 1 + k] = 0.;
+
+ u[3 + 2 * (order - 1) + k] = 0.;
+ v[3 + 2 * (order - 1) + k] = order - 1 - k;
+ w[3 + 2 * (order - 1) + k] = 0.;
}
- if (order >2){
- int nbdoftemp = nbdoftriangle(order-3);
- nodepositionface0(order-3,&u[3+3*(order-1)], &v[3+3*(order-1)], &w[3+3*(order-1)]);
- for (int k=0; k< nbdoftemp; k++){
- u[3+k+3*(order-1)] = u[3+k+3*(order-1)]*(order-3)+1.;
- v[3+k+3*(order-1)] = v[3+k+3*(order-1)]*(order-3)+1.;
- w[3+k+3*(order-1)] = w[3+k+3*(order-1)]*(order-3);
- }
+ if (order > 2){
+ int nbdoftemp = nbdoftriangle(order - 3);
+ nodepositionface0(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
+ &w[3 + 3* (order - 1)]);
+ for (int k = 0; k < nbdoftemp; k++){
+ u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3);
+ }
}
- for (int k=0; k< ndofT; k++){
- u[k] = u[k]/order;
- v[k] = v[k]/order;
- w[k] = w[k]/order;
+ for (int k = 0; k < ndofT; k++){
+ u[k] = u[k] / order;
+ v[k] = v[k] / order;
+ w[k] = w[k] / order;
}
- return;
-};
-
+}
-// -----------------------------------------------------------------------------
-void nodepositionface1(int order, double *u, double *v, double *w){
+void nodepositionface1(int order, double *u, double *v, double *w)
+{
int ndofT = nbdoftriangle(order);
- if (order ==0) {u[0]=0.;v[0]=0.; return;}
+ if (order == 0) { u[0] = 0.; v[0] = 0.; return; }
u[0]= 0.; v[0]= 0.; w[0] = 0.;
u[1]= order; v[1]= 0; w[1] = 0.;
u[2]= 0.; v[2]= 0.; w[2] = order;
// edges
- for (int k=0; k< (order-1); k++){
- u[3+k] = k+1; v[3+k]=0.; w[3+k] = 0.;
- u[3+order-1+k] = order -1 -k; v[3+order-1+k] = 0.; w[3+order-1+k] = k+1;
- u[3+2*(order-1)+k] = 0. ; v[3+2*(order-1)+k] = 0.; w[3+2*(order-1)+k] = order -1 -k;
+ for (int k = 0; k < (order - 1); k++){
+ u[3 + k] = k + 1;
+ v[3 + k] = 0.;
+ w[3 + k] = 0.;
+
+ u[3 + order - 1 + k] = order - 1 - k;
+ v[3 + order - 1 + k] = 0.;
+ w[3 + order - 1+ k ] = k + 1;
+
+ u[3 + 2 * (order - 1) + k] = 0. ;
+ v[3 + 2 * (order - 1) + k] = 0.;
+ w[3 + 2 * (order - 1) + k] = order - 1 - k;
}
- if (order >2){
- int nbdoftemp = nbdoftriangle(order-3);
- nodepositionface1(order-3,&u[3+3*(order-1)], &v[3+3*(order-1)], &w[3+3*(order-1)]);
- for (int k=0; k< nbdoftemp; k++){
- u[3+k+3*(order-1)] = u[3+k+3*(order-1)]*(order-3)+1.;
- v[3+k+3*(order-1)] = v[3+k+3*(order-1)]*(order-3);
- w[3+k+3*(order-1)] = w[3+k+3*(order-1)]*(order-3)+1.;
- }
+ if (order > 2){
+ int nbdoftemp = nbdoftriangle(order - 3);
+ nodepositionface1(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order -1 )],
+ &w[3 + 3 * (order - 1)]);
+ for (int k = 0; k < nbdoftemp; k++){
+ u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3);
+ w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ }
}
- for (int k=0; k< ndofT; k++){
- u[k] = u[k]/order;
- v[k] = v[k]/order;
- w[k] = w[k]/order;
+ for (int k = 0; k < ndofT; k++){
+ u[k] = u[k] / order;
+ v[k] = v[k] / order;
+ w[k] = w[k] / order;
}
- return;
-};
+}
-// -----------------------------------------------------------------------------
-void nodepositionface2(int order, double *u, double *v, double *w){
+void nodepositionface2(int order, double *u, double *v, double *w)
+{
int ndofT = nbdoftriangle(order);
- if (order ==0) {u[0]=0.;v[0]=0.; return;}
+ if (order == 0) { u[0] = 0.; v[0] = 0.; return; }
- u[0]= order; v[0]= 0.; w[0] = 0.;
- u[1]= 0.; v[1]= order; w[1] = 0.;
- u[2]= 0.; v[2]= 0.; w[2] = order;
+ u[0]= order; v[0]= 0.; w[0] = 0.;
+ u[1]= 0.; v[1]= order; w[1] = 0.;
+ u[2]= 0.; v[2]= 0.; w[2] = order;
// edges
- for (int k=0; k< (order-1); k++){
- u[3+k] = order -1 -k; v[3+k]= k+1; w[3+k] = 0.;
- u[3+order-1+k] = 0.; v[3+order-1+k] = order -1 -k; w[3+order-1+k] = k+1;
- u[3+2*(order-1)+k] = k+1; v[3+2*(order-1)+k] = 0.; w[3+2*(order-1)+k] = order -1 -k;
+ for (int k = 0; k < (order - 1); k++){
+ u[3 + k] = order - 1 - k;
+ v[3 + k] = k + 1;
+ w[3 + k] = 0.;
+
+ u[3 + order - 1 + k] = 0.;
+ v[3 + order - 1 + k] = order - 1 - k;
+ w[3 + order - 1 + k] = k + 1;
+
+ u[3 + 2 * (order - 1) + k] = k + 1;
+ v[3 + 2 * (order - 1) + k] = 0.;
+ w[3 + 2 * (order - 1) + k] = order - 1 - k;
}
- if (order >2){
- int nbdoftemp = nbdoftriangle(order-3);
- nodepositionface2(order-3,&u[3+3*(order-1)], &v[3+3*(order-1)], &w[3+3*(order-1)]);
- for (int k=0; k< nbdoftemp; k++){
- u[3+k+3*(order-1)] = u[3+k+3*(order-1)]*(order-3)+1.;
- v[3+k+3*(order-1)] = v[3+k+3*(order-1)]*(order-3)+1.;
- w[3+k+3*(order-1)] = w[3+k+3*(order-1)]*(order-3)+1.;
- }
+ if (order > 2){
+ int nbdoftemp = nbdoftriangle(order - 3);
+ nodepositionface2(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
+ &w[3 + 3 * (order - 1)]);
+ for (int k = 0; k < nbdoftemp; k++){
+ u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ }
}
- for (int k=0; k< ndofT; k++){
- u[k] = u[k]/order;
- v[k] = v[k]/order;
- w[k] = w[k]/order;
+ for (int k = 0; k < ndofT; k++){
+ u[k] = u[k] / order;
+ v[k] = v[k] / order;
+ w[k] = w[k] / order;
}
- return;
-};
+}
-// -----------------------------------------------------------------------------
-void nodepositionface3(int order, double *u, double *v, double *w){
+void nodepositionface3(int order, double *u, double *v, double *w)
+{
int ndofT = nbdoftriangle(order);
- if (order ==0) {u[0]=0.;v[0]=0.; return;}
+ if (order == 0) { u[0] = 0.; v[0] = 0.; return; }
- u[0]= 0.; v[0]= 0.; w[0] = 0.;
- u[1]= 0.; v[1]= order; w[1] = 0.;
- u[2]= 0.; v[2]= 0.; w[2] = order;
+ u[0]= 0.; v[0]= 0.; w[0] = 0.;
+ u[1]= 0.; v[1]= order; w[1] = 0.;
+ u[2]= 0.; v[2]= 0.; w[2] = order;
// edges
- for (int k=0; k< (order-1); k++){
- u[3+k] = 0.; v[3+k]= k+1; w[3+k] = 0.;
- u[3+order-1+k] = 0.; v[3+order-1+k] = order -1 -k; w[3+order-1+k] = k+1;
- u[3+2*(order-1)+k] = 0.; v[3+2*(order-1)+k] = 0.; w[3+2*(order-1)+k] = order -1 -k;
+ for (int k = 0; k < (order - 1); k++){
+ u[3 + k] = 0.;
+ v[3 + k]= k + 1;
+ w[3 + k] = 0.;
+
+ u[3 + order - 1 + k] = 0.;
+ v[3 + order - 1 + k] = order - 1 - k;
+ w[3 + order - 1 + k] = k + 1;
+
+ u[3 + 2 * (order - 1) + k] = 0.;
+ v[3 + 2 * (order - 1) + k] = 0.;
+ w[3 + 2 * (order - 1) + k] = order - 1 - k;
}
- if (order >2){
- int nbdoftemp = nbdoftriangle(order-3);
- nodepositionface3(order-3,&u[3+3*(order-1)], &v[3+3*(order-1)], &w[3+3*(order-1)]);
- for (int k=0; k< nbdoftemp; k++){
- u[3+k+3*(order-1)] = u[3+k+3*(order-1)]*(order-3);
- v[3+k+3*(order-1)] = v[3+k+3*(order-1)]*(order-3)+1.;
- w[3+k+3*(order-1)] = w[3+k+3*(order-1)]*(order-3)+1.;
- }
+ if (order > 2){
+ int nbdoftemp = nbdoftriangle(order - 3);
+ nodepositionface3(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
+ &w[3 + 3 * (order - 1)]);
+ for (int k = 0; k < nbdoftemp; k++){
+ u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3);
+ v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ }
}
- for (int k=0; k< ndofT; k++){
- u[k] = u[k]/order;
- v[k] = v[k]/order;
- w[k] = w[k]/order;
+ for (int k = 0; k < ndofT; k++){
+ u[k] = u[k] / order;
+ v[k] = v[k] / order;
+ w[k] = w[k] / order;
}
- return;
}
-// -----------------------------------------------------------------------------
-Double_Matrix gmshGeneratePointsTetrahedron(int order,bool serendip) {
-
- int nbPoints = (serendip ?
- 4 + 6 * std::max(0,order-1) + 4 * std::max(0,(order-2)*(order-1)/2):
- (order+1)*(order+2)*(order+3)/6);
+Double_Matrix gmshGeneratePointsTetrahedron(int order, bool serendip)
+{
+ int nbPoints =
+ (serendip ?
+ 4 + 6 * std::max(0, order - 1) + 4 * std::max(0, (order - 2) * (order - 1) / 2) :
+ (order + 1) * (order + 2) * (order + 3) / 6);
- Double_Matrix point(nbPoints,3);
+ Double_Matrix point(nbPoints, 3);
- double overOrder = 1./order;
+ double overOrder = 1. / order;
- point(0,0)= 0.;
- point(0,1)= 0.;
- point(0,2)= 0.;
+ point(0, 0) = 0.;
+ point(0, 1) = 0.;
+ point(0, 2) = 0.;
if (order > 0) {
+ point(1, 0) = order;
+ point(1, 1) = 0;
+ point(1, 2) = 0;
-
- point(1,0)= order;
- point(1,1)= 0;
- point(1,2)= 0;
+ point(2, 0) = 0.;
+ point(2, 1) = order;
+ point(2, 2) = 0.;
- point(2,0)= 0.;
- point(2,1)= order;
- point(2,2)= 0.;
-
- point(3,0)= 0.;
- point(3,1)= 0.;
- point(3,2)= order;
+ point(3, 0) = 0.;
+ point(3, 1) = 0.;
+ point(3, 2) = order;
if (order > 1) {
- for (int k=0; k< (order-1); k++) {
- point(4+k,0) = k+1;
- point(4+order-1+k,0) = order-1-k;
- point(4+2*(order-1)+k,0) = 0.;
- point(4+3*(order-1)+k,0) = 0.;
- point(4+4*(order-1)+k,0) = order-1-k;
- point(4+5*(order-1)+k,0) = 0.;
+ for (int k = 0; k < (order - 1); k++) {
+ point(4 + k, 0) = k + 1;
+ point(4 + order - 1 + k, 0) = order - 1 - k;
+ point(4 + 2 * (order - 1) + k, 0) = 0.;
+ point(4 + 3 * (order - 1) + k, 0) = 0.;
+ point(4 + 4 * (order - 1) + k, 0) = order - 1 - k;
+ point(4 + 5 * (order - 1) + k, 0) = 0.;
- point(4+k,1) = 0.;
- point(4+order-1+k,1) = k+1;
- point(4+2*(order-1)+k,1) = order-1-k;
- point(4+3*(order-1)+k,1) = 0.;
- point(4+4*(order-1)+k,1) = 0.;
- point(4+5*(order-1)+k,1) = order-1-k;
+ point(4 + k, 1) = 0.;
+ point(4 + order - 1 + k, 1) = k + 1;
+ point(4 + 2 * (order - 1) + k, 1) = order - 1 - k;
+ point(4 + 3 * (order - 1) + k, 1) = 0.;
+ point(4 + 4 * (order - 1) + k, 1) = 0.;
+ point(4 + 5 * (order - 1) + k, 1) = order - 1 - k;
- point(4+k,2) = 0.;
- point(4+order-1+k,2) = 0.;
- point(4+2*(order-1)+k,2) = 0.;
- point(4+3*(order-1)+k,2) = k+1;
- point(4+4*(order-1)+k,2) = k+1;
- point(4+5*(order-1)+k,2) = k+1;
+ point(4 + k, 2) = 0.;
+ point(4 + order - 1 + k, 2) = 0.;
+ point(4 + 2 * (order - 1) + k, 2) = 0.;
+ point(4 + 3 * (order - 1) + k, 2) = k + 1;
+ point(4 + 4 * (order - 1) + k, 2) = k + 1;
+ point(4 + 5 * (order - 1) + k, 2) = k + 1;
}
if (order > 2) {
- int ns = 4+6*(order-1);
- int nbdofface = nbdoftriangle(order-3);
+ int ns = 4 + 6 * (order - 1);
+ int nbdofface = nbdoftriangle(order - 3);
double *u = new double[nbdofface];
double *v = new double[nbdofface];
double *w = new double[nbdofface];
- nodepositionface0(order-3, u,v,w);
+ nodepositionface0(order - 3, u, v, w);
- for (int i=0; i < nbdofface; i++){
- point(ns+i,0) = u[i]*(order-3) + 1.;
- point(ns+i,1) = v[i]*(order-3) + 1.;
- point(ns+i,2) = w[i]*(order-3);
+ for (int i = 0; i < nbdofface; i++){
+ point(ns + i, 0) = u[i] * (order - 3) + 1.;
+ point(ns + i, 1) = v[i] * (order - 3) + 1.;
+ point(ns + i, 2) = w[i] * (order - 3);
}
ns = ns + nbdofface;
- nodepositionface1(order-3, u,v,w);
+ nodepositionface1(order - 3, u, v, w);
for (int i=0; i < nbdofface; i++){
- point(ns+i,0) = u[i]*(order-3) + 1.;
- point(ns+i,1) = v[i]*(order-3) ;
- point(ns+i,2) = w[i]*(order-3) + 1.;
+ point(ns + i, 0) = u[i] * (order - 3) + 1.;
+ point(ns + i, 1) = v[i] * (order - 3) ;
+ point(ns + i, 2) = w[i] * (order - 3) + 1.;
}
ns = ns + nbdofface;
- nodepositionface2(order-3, u,v,w);
+ nodepositionface2(order - 3, u, v, w);
- for (int i=0; i < nbdofface; i++){
- point(ns+i,0) = u[i]*(order-3) + 1.;
- point(ns+i,1) = v[i]*(order-3) + 1.;
- point(ns+i,2) = w[i]*(order-3) + 1.;
+ for (int i = 0; i < nbdofface; i++){
+ point(ns + i, 0) = u[i] * (order - 3) + 1.;
+ point(ns + i, 1) = v[i] * (order - 3) + 1.;
+ point(ns + i, 2) = w[i] * (order - 3) + 1.;
}
- ns =ns +nbdofface;
+ ns = ns + nbdofface;
- nodepositionface3(order-3, u,v,w);
+ nodepositionface3(order - 3, u, v, w);
- for (int i=0; i < nbdofface; i++){
- point(ns+i,0) = u[i]*(order-3);
- point(ns+i,1) = v[i]*(order-3) + 1.;
- point(ns+i,2) = w[i]*(order-3) + 1.;
+ for (int i = 0; i < nbdofface; i++){
+ point(ns + i, 0) = u[i] * (order - 3);
+ point(ns + i, 1) = v[i] * (order - 3) + 1.;
+ point(ns + i, 2) = w[i] * (order - 3) + 1.;
}
- ns =ns + nbdofface;
+ ns = ns + nbdofface;
delete [] u;
delete [] v;
delete [] w;
-
- if (!serendip && order > 3) {
+ if (!serendip && order > 3) {
- Double_Matrix interior = gmshGeneratePointsTetrahedron (order-4,false); // tetrahedron order - 4
- for (size_t k=0; k< interior.size1() ; k++) {
- point(ns+k,0) = 1.+interior(k,0)*(order-4);
- point(ns+k,1) = 1.+interior(k,1)*(order-4);
- point(ns+k,2) = 1.+interior(k,2)*(order-4);
+ Double_Matrix interior = gmshGeneratePointsTetrahedron(order - 4, false);
+ for (size_t k = 0; k < interior.size1(); k++) {
+ point(ns + k, 0) = 1. + interior(k, 0) * (order - 4);
+ point(ns + k, 1) = 1. + interior(k, 1) * (order - 4);
+ point(ns + k, 2) = 1. + interior(k, 2) * (order - 4);
}
}
}
@@ -400,23 +437,21 @@ Double_Matrix gmshGeneratePointsTetrahedron(int order,bool serendip) {
return point;
}
-
-Double_Matrix gmshGeneratePointsTriangle(int order,bool serendip) {
-
- int nbPoints = serendip ? 3 * order : (order+1)*(order+2)/2;
- Double_Matrix point(nbPoints,2);
+Double_Matrix gmshGeneratePointsTriangle(int order, bool serendip)
+{
+ int nbPoints = serendip ? 3 * order : (order + 1) * (order + 2) / 2;
+ Double_Matrix point(nbPoints, 2);
- point(0,0) = 0;
- point(0,1) = 0;
+ point(0, 0) = 0;
+ point(0, 1) = 0;
- double dd = 1./order;
+ double dd = 1. / order;
if (order > 0) {
-
- point(1,0) = 1;
- point(1,1) = 0;
- point(2,0) = 0;
- point(2,1) = 1;
+ point(1, 0) = 1;
+ point(1, 1) = 0;
+ point(2, 0) = 0;
+ point(2, 1) = 1;
int index = 3;
@@ -425,56 +460,57 @@ Double_Matrix gmshGeneratePointsTriangle(int order,bool serendip) {
double ksi = 0;
double eta = 0;
- for (int i=0;i<order-1;i++,index++) {
+ for (int i = 0; i < order - 1; i++, index++) {
ksi += dd;
- point(index,0) = ksi;
- point(index,1) = eta;
+ point(index, 0) = ksi;
+ point(index, 1) = eta;
}
ksi = 1.;
- for (int i=0;i<order-1;i++,index++) {
+ for (int i = 0; i < order - 1; i++, index++) {
ksi -= dd;
eta += dd;
- point(index,0) = ksi;
- point(index,1) = eta;
+ point(index, 0) = ksi;
+ point(index, 1) = eta;
}
eta = 1.;
ksi = 0.;
- for (int i=0;i<order-1;i++,index++) {
+ for (int i = 0; i < order - 1; i++, index++) {
eta -= dd;
- point(index,0) = ksi;
- point(index,1) = eta;
+ point(index, 0) = ksi;
+ point(index, 1) = eta;
}
if (order > 2 && !serendip) {
Double_Matrix inner = gmshGeneratePointsTriangle(order - 3, serendip);
- inner.scale(1.-3.*dd);
+ inner.scale(1. - 3. * dd);
inner.add(dd);
- Double_Matrix(point.touchSubmatrix(index,nbPoints - index,0,2)).memcpy(inner);
+ Double_Matrix(point.touchSubmatrix(index, nbPoints - index, 0, 2)).memcpy(inner);
}
}
}
return point;
}
-Double_Matrix generateLagrangeMonomialCoefficients(const Double_Matrix& monomial,const Double_Matrix& point) {
-
+Double_Matrix generateLagrangeMonomialCoefficients(const Double_Matrix& monomial,
+ const Double_Matrix& point)
+{
if (monomial.size1() != point.size1()) throw;
if (monomial.size2() != point.size2()) throw;
int ndofs = monomial.size1();
int dim = monomial.size2();
- Double_Matrix Vandermonde(ndofs,ndofs);
+ Double_Matrix Vandermonde(ndofs, ndofs);
- for (int i=0;i<ndofs;i++) {
- for (int j=0;j<ndofs;j++) {
+ for (int i = 0; i < ndofs; i++) {
+ for (int j = 0; j < ndofs; j++) {
double dd = 1.;
- for (int k=0;k<dim;k++) dd *= pow(point(i,k),monomial(j,k));
- Vandermonde(i,j) = dd;
+ for (int k = 0; k < dim; k++) dd *= pow(point(i, k), monomial(j, k));
+ Vandermonde(i, j) = dd;
}
}
@@ -482,78 +518,78 @@ Double_Matrix generateLagrangeMonomialCoefficients(const Double_Matrix& monomial
double det = Vandermonde.determinant();
- if (det == 0.0)throw;
+ if (det == 0.0) throw;
- Double_Matrix coefficient(ndofs,ndofs);
+ Double_Matrix coefficient(ndofs, ndofs);
- for (int i=0;i<ndofs;i++) {
- for (int j=0;j<ndofs;j++) {
- int f = (i+j)%2 == 0 ? 1 : -1;
- Double_Matrix cofactor = Vandermonde.cofactor(i,j);
- coefficient(i,j) = f * cofactor.determinant()/det;
+ for (int i = 0; i < ndofs; i++) {
+ for (int j = 0; j < ndofs; j++) {
+ int f = (i + j) % 2 == 0 ? 1 : -1;
+ Double_Matrix cofactor = Vandermonde.cofactor(i, j);
+ coefficient(i, j) = f * cofactor.determinant() / det;
}
}
Vandermonde.set_all(0.);
- for (int i=0;i<ndofs;i++) {
- for (int j=0;j<ndofs;j++) {
+ for (int i = 0; i < ndofs; i++) {
+ for (int j = 0; j < ndofs; j++) {
double dd = 1.;
- for (int k=0;k<dim;k++) dd *= pow(point(i,k),monomial(j,k));
-
- for (int k=0;k<ndofs;k++) {
- Vandermonde(i,k) += coefficient(k,j)*dd;
+ for (int k = 0; k < dim; k++) dd *= pow(point(i, k), monomial(j, k));
+ for (int k = 0; k < ndofs; k++) {
+ Vandermonde(i, k) += coefficient(k, j) * dd;
}
}
}
return coefficient;
}
-std::map<int,gmshFunctionSpace> gmshFunctionSpaces::fs ;
+std::map<int, gmshFunctionSpace> gmshFunctionSpaces::fs;
-const gmshFunctionSpace & gmshFunctionSpaces::find (int tag) {
+const gmshFunctionSpace &gmshFunctionSpaces::find(int tag)
+{
std::map<int,gmshFunctionSpace>::const_iterator it = fs.find(tag);
- if (it != fs.end())return it->second;
+ if (it != fs.end()) return it->second;
gmshFunctionSpace F;
- switch ( tag ){
+ switch (tag){
case MSH_TRI_3 :
- F.monomials = generatePascalTriangle (1);
+ F.monomials = generatePascalTriangle(1);
F.points = gmshGeneratePointsTriangle(1, false);
break;
case MSH_TRI_6 :
- F.monomials = generatePascalTriangle (2);
+ F.monomials = generatePascalTriangle(2);
F.points = gmshGeneratePointsTriangle(2, false);
break;
case MSH_TRI_9 :
- F.monomials = generatePascalSerendipityTriangle (3);
+ F.monomials = generatePascalSerendipityTriangle(3);
F.points = gmshGeneratePointsTriangle(3, true);
break;
case MSH_TRI_10 :
- F.monomials = generatePascalTriangle (3);
+ F.monomials = generatePascalTriangle(3);
F.points = gmshGeneratePointsTriangle(3, false);
break;
case MSH_TRI_12 :
- F.monomials = generatePascalSerendipityTriangle (4);
+ F.monomials = generatePascalSerendipityTriangle(4);
F.points = gmshGeneratePointsTriangle(4, true);
break;
case MSH_TRI_15 :
- F.monomials = generatePascalTriangle (4);
+ F.monomials = generatePascalTriangle(4);
F.points = gmshGeneratePointsTriangle(4, false);
break;
case MSH_TRI_15I :
- F.monomials = generatePascalSerendipityTriangle (5);
+ F.monomials = generatePascalSerendipityTriangle(5);
F.points = gmshGeneratePointsTriangle(5, true);
break;
case MSH_TRI_21 :
- F.monomials = generatePascalTriangle (5);
+ F.monomials = generatePascalTriangle(5);
F.points = gmshGeneratePointsTriangle(5, false);
break;
default :
throw;
}
- F.coefficients = generateLagrangeMonomialCoefficients ( F.monomials, F.points );
- fs.insert(std::make_pair(tag,F));
+ F.coefficients = generateLagrangeMonomialCoefficients(F.monomials, F.points);
+ fs.insert(std::make_pair(tag, F));
return fs[tag];
}
diff --git a/Numeric/FunctionSpace.h b/Numeric/FunctionSpace.h
index aebad28901b8226bce555779829bc70841ae79d7..308313207123c795ecae3a76dfec90e4d4e755ff 100644
--- a/Numeric/FunctionSpace.h
+++ b/Numeric/FunctionSpace.h
@@ -1,7 +1,7 @@
-#ifndef _GMSH_FCTSP_H_
-#define _GMSH_FCTSP_H_
+#ifndef _FUNCTION_SPACE_H_
+#define _FUNCTION_SPACE_H_
-// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
+// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
@@ -19,80 +19,80 @@
// USA.
//
// Please report all bugs and problems to <gmsh@geuz.org>.
+
#include <map>
#include "GmshMatrix.h"
+
struct gmshFunctionSpace
{
Double_Matrix points;
Double_Matrix monomials;
Double_Matrix coefficients;
- inline void computePows ( double uu, double vv , double p[][2]) const{
+ inline void computePows(double uu, double vv, double p[][2]) const
+ {
for(int j = 0; j < coefficients.size2(); j++){
- p[j][0] = pow(uu,monomials(j,0));
- p[j][1] = pow(vv,monomials(j,1));
+ p[j][0] = pow(uu,monomials(j, 0));
+ p[j][1] = pow(vv,monomials(j, 1));
}
}
- inline void f ( double u, double v, double w, double *sf) const{
- for (int i = 0; i < coefficients.size1(); i++){
+ inline void f(double u, double v, double w, double *sf) const
+ {
+ for(int i = 0; i < coefficients.size1(); i++){
sf[i] = 0;
for(int j = 0; j < coefficients.size2(); j++){
- sf[i] += coefficients(i,j) *
- pow(u,monomials(j,0)) *
- pow(v,monomials(j,1)) *
- pow(w,monomials(j,2));
+ sf[i] += coefficients(i, j) * pow(u, monomials(j, 0)) *
+ pow(v, monomials(j, 1)) * pow(w, monomials(j, 2));
}
}
}
- inline void f ( double u, double v, double *sf) const{
+ inline void f(double u, double v, double *sf) const
+ {
double p[256][2];
- computePows ( u, v , p);
+ computePows(u, v, p);
for (int i = 0; i < coefficients.size1(); i++){
sf[i] = 0;
for(int j = 0; j < coefficients.size2(); j++){
- sf[i] += coefficients(i,j) * p[j][0] * p[j][1] ;
+ sf[i] += coefficients(i, j) * p[j][0] * p[j][1];
}
}
}
- inline void df ( double u, double v, double w, double grads[][3]) const{
+ inline void df(double u, double v, double w, double grads[][3]) const
+ {
for (int i = 0; i < coefficients.size1(); i++){
grads[i][0] = 0;
grads[i][1] = 0;
grads[i][2] = 0;
for(int j = 0; j < coefficients.size2(); j++){
- if ((monomials)(j,0) > 0)
- grads[i][0] += (coefficients)(i,j) *
- pow(u,(monomials)(j,0)-1) * (monomials)(j,0) *
- pow(v,(monomials)(j,1)) *
- pow(w,(monomials)(j,2));
- if ((monomials)(j,1) > 0)
- grads[i][1] += (coefficients)(i,j) *
- pow(u,(monomials)(j,0)) *
- pow(v,(monomials)(j,1)-1) * (monomials)(j,1) *
- pow(w,(monomials)(j,2));
- if ((monomials)(j,2) > 0)
- grads[i][2] += (coefficients)(i,j) *
- pow(u,(monomials)(j,0)) *
- pow(v,(monomials)(j,1)) *
- pow(w,(monomials)(j,2)-1) * (monomials)(j,2) ;
+ if ((monomials)(j, 0) > 0)
+ grads[i][0] += (coefficients)(i, j) * pow(u, (monomials)(j, 0) - 1) *
+ (monomials)(j, 0) * pow(v, (monomials)(j, 1)) *
+ pow(w, (monomials)(j, 2));
+ if ((monomials)(j, 1) > 0)
+ grads[i][1] += (coefficients)(i, j) * pow(u,(monomials)(j, 0)) *
+ pow(v, (monomials)(j, 1) - 1) * (monomials)(j, 1) *
+ pow(w, (monomials)(j, 2));
+ if ((monomials)(j, 2) > 0)
+ grads[i][2] += (coefficients)(i, j) * pow(u, (monomials)(j, 0)) *
+ pow(v, (monomials)(j, 1)) * pow(w, (monomials)(j, 2) - 1) *
+ (monomials)(j, 2);
}
}
}
- inline void df ( double u, double v, double grads[][2]) const{
-
+ inline void df(double u, double v, double grads[][2]) const
+ {
double p[256][2];
- computePows ( u, v , p);
-
+ computePows(u, v, p);
for (int i = 0; i < coefficients.size1(); i++){
grads[i][0] = 0;
grads[i][1] = 0;
for(int j = 0; j < coefficients.size2(); j++){
- if ((monomials)(j,0) > 0)
- grads[i][0] += (coefficients)(i,j) *
- pow(u,(monomials)(j,0)-1) * (monomials)(j,0) * p[j][1];
- if ((monomials)(j,1) > 0)
- grads[i][1] += (coefficients)(i,j) * p[j][0] *
- pow(v,(monomials)(j,1)-1) * (monomials)(j,1);
+ if ((monomials)(j, 0) > 0)
+ grads[i][0] += (coefficients)(i, j) *
+ pow(u, (monomials)(j, 0) - 1) * (monomials)(j, 0) * p[j][1];
+ if ((monomials)(j, 1) > 0)
+ grads[i][1] += (coefficients)(i, j) * p[j][0] *
+ pow(v, (monomials)(j, 1) - 1) * (monomials)(j, 1);
}
}
}
@@ -100,10 +100,9 @@ struct gmshFunctionSpace
class gmshFunctionSpaces
{
- static std::map<int,gmshFunctionSpace> fs;
+ static std::map<int, gmshFunctionSpace> fs;
public :
- static const gmshFunctionSpace & find ( int );
+ static const gmshFunctionSpace &find(int);
};
-
#endif
diff --git a/Numeric/GaussLegendre1D.h b/Numeric/GaussLegendre1D.h
index 3d1de525479a7c893c0be9f74bd60bbdd05cd76d..a3da4ddd3e6f5c7cff6e957b07e9dcafc0f371c2 100644
--- a/Numeric/GaussLegendre1D.h
+++ b/Numeric/GaussLegendre1D.h
@@ -1,6 +1,7 @@
-#ifndef _GMSH_GL1D_
-#define _GMSH_GL1D_
-// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
+#ifndef _GMSH_GAUSS_LEGENDRE_1D_H_
+#define _GMSH_GAUSS_LEGENDRE_1D_H_
+
+// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
@@ -19,55 +20,56 @@
//
// Please report all bugs and problems to <gmsh@geuz.org>.
- /* 1 point rule points */
- static double _GL_pt1[1]={
+/* 1 point rule points */
+static double _GL_pt1[1]={
0.000000000000000e+00};
- /* 1 point rule weights */
- static double _GL_wt1[1]={
+/* 1 point rule weights */
+static double _GL_wt1[1]={
2.000000000000000e+00};
- /* 2 point rule points */
- static double _GL_pt2[2]={
- -5.773502691896257e-01, 5.773502691896257e-01};
+/* 2 point rule points */
+static double _GL_pt2[2]={
+ -5.773502691896257e-01, 5.773502691896257e-01};
- /* 2 point rule weights */
- static double _GL_wt2[2]={
+/* 2 point rule weights */
+static double _GL_wt2[2]={
1.000000000000000e+00, 1.000000000000000e+00};
- /* 3 point rule points */
- static double _GL_pt3[3]={
- -7.745966692414834e-01, 0.000000000000000e+00, 7.745966692414834e-01};
+/* 3 point rule points */
+static double _GL_pt3[3]={
+ -7.745966692414834e-01, 0.000000000000000e+00, 7.745966692414834e-01};
- /* 3 point rule weights */
- static double _GL_wt3[3]={
+/* 3 point rule weights */
+static double _GL_wt3[3]={
5.555555555555552e-01, 8.888888888888888e-01, 5.555555555555552e-01};
- /* 4 point rule points */
- static double _GL_pt4[4]={
- -8.611363115940526e-01,-3.399810435848563e-01, 3.399810435848563e-01, 8.611363115940526e-01};
+/* 4 point rule points */
+static double _GL_pt4[4]={
+ -8.611363115940526e-01,-3.399810435848563e-01, 3.399810435848563e-01, 8.611363115940526e-01};
- /* 4 point rule weights */
- static double _GL_wt4[4]={
+/* 4 point rule weights */
+static double _GL_wt4[4]={
3.478548451374537e-01, 6.521451548625464e-01, 6.521451548625464e-01, 3.478548451374537e-01};
- /* 5 point rule points */
- static double _GL_pt5[5]={
- -9.061798459386640e-01,-5.384693101056831e-01, 0.000000000000000e+00, 5.384693101056831e-01, 9.061798459386640e-01};
+/* 5 point rule points */
+static double _GL_pt5[5]={
+ -9.061798459386640e-01,-5.384693101056831e-01, 0.000000000000000e+00, 5.384693101056831e-01, 9.061798459386640e-01};
- /* 5 point rule weights */
- static double _GL_wt5[5]={
+/* 5 point rule weights */
+static double _GL_wt5[5]={
2.369268850561890e-01, 4.786286704993665e-01, 5.688888888888889e-01, 4.786286704993665e-01, 2.369268850561890e-01};
- /* 6 point rule points */
- static double _GL_pt6[6]={
- -9.324695142031521e-01,-6.612093864662646e-01,-2.386191860831969e-01, 2.386191860831969e-01, 6.612093864662646e-01, 9.324695142031521e-01};
+/* 6 point rule points */
+static double _GL_pt6[6]={
+ -9.324695142031521e-01,-6.612093864662646e-01,-2.386191860831969e-01, 2.386191860831969e-01, 6.612093864662646e-01, 9.324695142031521e-01};
- /* 6 point rule weights */
- static double _GL_wt6[6]={
+/* 6 point rule weights */
+static double _GL_wt6[6]={
1.713244923791705e-01, 3.607615730481386e-01, 4.679139345726913e-01, 4.679139345726913e-01, 3.607615730481386e-01, 1.713244923791705e-01};
-inline void gmshGaussLegendre1D (int nbQuadPoints , double **t, double **w){
+inline void gmshGaussLegendre1D(int nbQuadPoints , double **t, double **w)
+{
switch (nbQuadPoints){
case 1:
*t = _GL_pt1;