diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index ca677e92320a707e95a30971da652a06d223b4e0..d53fa4256420e9185af152b97568ac0fec5ed267 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -118,6 +118,9 @@ class GEdge : public GEntity {
// Resets the mesh attributes to default values
virtual void resetMeshAttributes();
+ // True if entity is periodic in the "dim" direction.
+ virtual bool periodic(int dim) const { return v0 == v1 ; }
+
struct {
char Method;
double coeffTransfinite;
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 8bb04cdcf9a5b93e189651c23351ce11090a2c8c..28dc03b1a7f0b27b18a3ad78606bcf3f4cc0cd41 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -1,4 +1,4 @@
-// $Id: BackgroundMesh.cpp,v 1.38 2008-02-21 12:11:12 geuzaine Exp $
+// $Id: BackgroundMesh.cpp,v 1.39 2008-03-12 08:36:48 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -118,9 +118,9 @@ double LC_MVertex_CURV(GEntity *ge, double U, double V)
case 1:
{
GEdge *ged = (GEdge *)ge;
- // Crv = ged->curvature(U);
- // Crv = std::max(Crv, max_surf_curvature(ged, U));
- 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;
case 2:
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 562e77dfbc89db2a2eb92ce6f75db821d7e7bace..ddec4fe8cd118f4ff60dd1ef99d88954c6ea77af 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1,4 +1,4 @@
-// $Id: HighOrder.cpp,v 1.24 2008-03-01 01:32:03 geuzaine Exp $
+// $Id: HighOrder.cpp,v 1.25 2008-03-12 08:36:49 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -78,7 +78,7 @@ bool reparamOnFace(MVertex *v, GFace *gf, SPoint2 ¶m)
static double mylength(GEdge *ge, int i, double *u)
{
- return ge->length(u[i], u[i+1], 6);
+ return ge->length(u[i], u[i+1], 10);
}
static void myresid(int N, GEdge *ge, double *u, Double_Vector &r)
@@ -88,7 +88,7 @@ static void myresid(int N, GEdge *ge, double *u, Double_Vector &r)
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)
+bool computeEquidistantParameters(GEdge *ge, double u0, double uN, int N, double *u, double underRelax)
{
const double PRECISION = 1.e-6;
const int MAX_ITER = 50;
@@ -98,8 +98,6 @@ bool computeEquidistantParameters(GEdge *ge, double u0, double uN, int N, double
// N is the total number of points (3 for quadratic, 4 for cubic ...)
// u[0] = u0;
// u[N-1] = uN;
- // dist(ge(u[i]), ge(u[i+1])) = dist(ge(u[i+1]), ge(u[i+2])), i = 0,...,N-2
-
// initialize as equidistant in parameter space
u[0] = u0;
double du = (uN - u0) / (N - 1);
@@ -135,24 +133,29 @@ bool computeEquidistantParameters(GEdge *ge, double u0, double uN, int N, double
J.lu_solve(R, DU);
for (int i = 0; i < M; i++){
- u[i+1] -= DU(i);
+ u[i+1] -= underRelax*DU(i);
}
+ // printf("N %d M %d u1 = %g u0 = %g uN1 = %22.15E uN = %22.15E\n",N,M,u[1],u0,u[N - 1],uN);
+
if (u[1] < u0) break;
- if (u[N - 1] > uN) break;
+ if (u[N - 2] > uN) break;
double newt_norm = DU.norm();
- if (newt_norm < PRECISION) return true;
+ // printf("%22.15E\n",newt_norm);
+ if (newt_norm < PRECISION) {/*printf("ok %g\n",underRelax);*/return true;}
}
// FAILED, use equidistant in param space
- for (int i = 1; i < N; i++){
- u[i] = u[i - 1] + du;
- }
+ // printf("coucou FAILED\n");
+ // printf("failed %g\n",underRelax);
+// 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);
+ return gf->length(SPoint2(u[i], v[i]), SPoint2(u[i + 1], v[i + 1]), 10);
}
static void myresid(int N, GFace *gf, double *u, double *v, Double_Vector &r)
@@ -267,17 +270,30 @@ void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
}
else{
MVertex *v0 = edge.getVertex(0), *v1 = edge.getVertex(1);
+
double u0 = 0., u1 = 0.;
bool reparamOK = true;
if(!linear && ge->geomType() != GEntity::DiscreteCurve &&
ge->geomType() != GEntity::BoundaryLayerCurve){
reparamOK &= reparamOnEdge(v0, ge, u0);
- reparamOK &= reparamOnEdge(v1, ge, u1);
+
+ if (ge->periodic(0) && v1 == ge->getEndVertex()->mesh_vertices[0]){
+ Range<double> par = ge->parBounds(0);
+ u1 = par.high();
+ }
+ else
+ reparamOK &= reparamOnEdge(v1, ge, u1);
}
double US[100];
if(reparamOK && !linear && ge->geomType() != GEntity::DiscreteCurve){
- computeEquidistantParameters(ge, u0, u1, nPts + 2, US);
- }
+ double relax = 1.;
+ while (1){
+ if (computeEquidistantParameters(ge, u0, u1, nPts + 2, US,relax))break;
+ relax /= 2.0;
+ if (relax < 1.e-2)break;
+ }
+ if (relax < 1.e-2)printf("failure in computing equidistant parameters %g\n",relax);
+ }
std::vector<MVertex*> temp;
for(int j = 0; j < nPts; j++){
const double t = (double)(j + 1)/(nPts + 1);
@@ -769,6 +785,32 @@ bool straightLine(std::vector<MVertex*> &l, MVertex *n1, MVertex *n2)
return true;
}
+static double mesh_functional_distorsion ( MTriangle *t, double u, double v){
+ // compute uncurved element jacobian d_u x and d_v x
+ double mat[2][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];
+ prodve(v1, v2, normal1);
+ double nn = sqrt(DSQR(normal1[0]) + DSQR(normal1[1]) + DSQR(normal1[2]));
+
+ // compute uncurved element jacobian d_u x and d_v x
+ t->jac(u, v, mat);
+ double v1b[3] = {mat[0][0], mat[0][1], mat[0][2]};
+ double v2b[3] = {mat[1][0], mat[1][1], mat[1][2]};
+ double normal[3];
+ prodve(v1b, v2b, normal);
+
+ double sign;
+ prosca(normal1, normal, &sign);
+ double det = norm3(normal) * (sign > 0 ? 1. : -1.) / nn;
+
+ // compute distorsion
+ double dist = std::min(1. / det, det);
+ return dist;
+}
+
void getMinMaxJac (MTriangle *t, double &minJ, double &maxJ)
{
double mat[2][3];
@@ -1186,9 +1228,9 @@ void checkHighOrderTriangles(GModel *m)
void printJacobians(GModel *m, const char *nm)
{
- return;
+ // return;
- const int n = 22;
+ const int n = 5;
double D[n][n];
double X[n][n];
double Y[n][n];
@@ -1199,17 +1241,13 @@ void printJacobians(GModel *m, const char *nm)
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
for(unsigned int i = 0; i < (*it)->triangles.size(); i++){
MTriangle *t = (*it)->triangles[i];
- double mat[2][3];
- t->jac(1,0,0,0,mat);
- double dj0 = mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0];
-
for(int i = 0; i < n; i++){
for(int k = 0; k < n - i; k++){
SPoint3 pt;
- t->jac((double)i / (n - 1), (double)k / (n - 1), mat);
- t->pnt((double)i / (n - 1), (double)k / (n - 1), pt);
- const double det = mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0];
- D[i][k] = std::min(det/dj0,dj0/det);
+ double u = (double)i / (n - 1);
+ double v = (double)k / (n - 1);
+ t->pnt(u,v, pt);
+ D[i][k] = mesh_functional_distorsion (t,u,v);
X[i][k] = pt.x();
Y[i][k] = pt.y();
Z[i][k] = pt.z();
@@ -1277,6 +1315,7 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
printJacobians(m, "detjIni.pos");
+
if(CTX.mesh.smooth_internal_edges){
checkHighOrderTriangles(m);
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
@@ -1296,6 +1335,7 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
}
printJacobians(m, "detjOpt.pos");
+
checkHighOrderTriangles(m);
double t2 = Cpu();
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 16438dfcfec42f4d63924f6822bbbf262e8bec30..0a3187b012ed35fee67df7f37e31c0968e2e717f 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.123 2008-03-01 01:32:03 geuzaine Exp $
+// $Id: meshGFace.cpp,v 1.124 2008-03-12 08:36:49 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -1275,8 +1275,15 @@ void deMeshGFace::operator() (GFace *gf)
gf->meshStatistics.nbTriangle = gf->meshStatistics.nbEdge = 0;
}
+const int debugSurface = -1;
+
void meshGFace::operator() (GFace *gf)
{
+ if (debugSurface>=0 && gf->tag() != debugSurface){
+ gf->meshStatistics.status = GFace::DONE;
+ return;
+ }
+
if(gf->geomType() == GEntity::DiscreteSurface) return;
if(gf->geomType() == GEntity::BoundaryLayerSurface) return;
if(gf->geomType() == GEntity::ProjectionFace) return;
@@ -1317,10 +1324,10 @@ void meshGFace::operator() (GFace *gf)
Msg(DEBUG1, "Generating the mesh");
if(noseam(gf) || gf->getNativeType() == GEntity::GmshModel || gf->edgeLoops.empty()){
- gmsh2DMeshGenerator(gf,0, false);
+ gmsh2DMeshGenerator(gf,0, debugSurface>=0);
}
else{
- if(!gmsh2DMeshGeneratorPeriodic(gf, false))
+ if(!gmsh2DMeshGeneratorPeriodic(gf, debugSurface>=0))
Msg(GERROR, "Impossible to mesh face %d", gf->tag());
}
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 657732f708247ed8dcbfdbe20e22b3278654c371..cca3ca113c041aa3f0bf865aa8ef23450971ad0e 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceDelaunayInsertion.cpp,v 1.11 2008-02-21 13:34:40 geuzaine Exp $
+// $Id: meshGFaceDelaunayInsertion.cpp,v 1.12 2008-03-12 08:36:49 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -467,7 +467,14 @@ void insertVerticesInFace(GFace *gf, BDS_Mesh *bds)
Vs[base->getVertex(2)->getNum()]) / 3.};
buildMetric(gf, pa, metric);
circumCenterMetric(worst->tri(), metric, Us, Vs, center, r2);
+
bool inside = invMapUV(worst->tri(), center, Us, Vs, uv, 1.e-8);
+ if (!inside && worst->getNeigh(0))
+ inside |= invMapUV(worst->getNeigh(0)->tri(), center, Us, Vs, uv, 1.e-8);
+ if (!inside && worst->getNeigh(1))
+ inside |= invMapUV(worst->getNeigh(1)->tri(), center, Us, Vs, uv, 1.e-8);
+ if (!inside && worst->getNeigh(2))
+ inside |= invMapUV(worst->getNeigh(2)->tri(), center, Us, Vs, uv, 1.e-8);
if (inside) {
// we use here local coordinates as real coordinates
// x,y and z will be computed hereafter
@@ -497,6 +504,13 @@ void insertVerticesInFace(GFace *gf, BDS_Mesh *bds)
gf->mesh_vertices.push_back(v);
}
else {
+ Msg(DEBUG,"Point %g %g is outside (%g %g , %g %g , %g %g) (metric %g %g %g)",center[0],center[1],
+ Us[base->getVertex(0)->getNum()],
+ Vs[base->getVertex(0)->getNum()],
+ Us[base->getVertex(1)->getNum()],
+ Vs[base->getVertex(1)->getNum()],
+ Us[base->getVertex(2)->getNum()],
+ Vs[base->getVertex(2)->getNum()],metric[0],metric[1],metric[2]);
AllTris.erase(AllTris.begin());
worst->forceRadius(0);
AllTris.insert(worst);
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 2c0e68dfcfcad2be0d01db3c8aa76e60d35f6c61..a6f8af177a4c1f936b3adff5787abe632fdc7ba3 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegionDelaunayInsertion.cpp,v 1.39 2008-03-06 14:19:01 geuzaine Exp $
+// $Id: meshGRegionDelaunayInsertion.cpp,v 1.40 2008-03-12 08:36:49 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -131,9 +131,12 @@ void recurFindCavity(std::list<faceXtet> & shell,
}
}
-bool insertVertex(MVertex *v, MTet4 *t, MTet4Factory &myFactory,
- std::set<MTet4*, compareTet4Ptr> &allTets,
- std::vector<double> &vSizes)
+bool insertVertex(MVertex *v,
+ MTet4 *t,
+ MTet4Factory &myFactory,
+ std::set<MTet4*,compareTet4Ptr> &allTets,
+ std::vector<double> & vSizes,
+ std::vector<double> & vSizesBGM)
{
std::list<faceXtet> shell;
std::list<MTet4*> cavity;
@@ -206,7 +209,7 @@ bool insertVertex(MVertex *v, MTet4 *t, MTet4Factory &myFactory,
// it->v[2]->y(),
// it->v[2]->z());
- MTet4 *t4 = myFactory.Create(tr, vSizes);
+ MTet4 *t4 = myFactory.Create(tr, vSizes, vSizesBGM);
t4->setOnWhat(t->onWhat());
newTets[k++] = t4;
// all new tets are pushed front in order to ba able to destroy
@@ -722,6 +725,7 @@ void insertVerticesInRegion (GRegion *gr)
// sizeof(MTet4), sizeof(MTetrahedron), sizeof(MVertex));
std::vector<double> vSizes;
+ std::vector<double> vSizesBGM;
MTet4Factory myFactory(1600000);
std::set<MTet4*, compareTet4Ptr> &allTets = myFactory.getAllTets();
int NUM = 0;
@@ -734,11 +738,12 @@ void insertVerticesInRegion (GRegion *gr)
it != vSizesMap.end(); ++it){
it->first->setNum(NUM++);
vSizes.push_back(it->second);
+ vSizesBGM.push_back(it->second);
}
}
for(unsigned int i = 0; i < gr->tetrahedra.size(); i++)
- allTets.insert(myFactory.Create(gr->tetrahedra[i], vSizes));
+ allTets.insert(myFactory.Create(gr->tetrahedra[i], vSizes,vSizesBGM));
gr->tetrahedra.clear();
connectTets(allTets.begin(), allTets.end());
@@ -815,10 +820,12 @@ void insertVerticesInRegion (GRegion *gr)
uvw[0] * vSizes[worst->tet()->getVertex(1)->getNum()] +
uvw[1] * vSizes[worst->tet()->getVertex(2)->getNum()] +
uvw[2] * vSizes[worst->tet()->getVertex(3)->getNum()];
- double lc = std::min(lc1, BGM_MeshSize(gr, 0, 0, center[0], center[1], center[2]));
- vSizes.push_back(lc);
+ double lc = BGM_MeshSize(gr, 0, 0, center[0], center[1], center[2]);
+ // double lc = std::min(lc1, BGM_MeshSize(gr, 0, 0, center[0], center[1], center[2]));
+ vSizes.push_back(lc1);
+ vSizesBGM.push_back(lc);
// compute mesh spacing there
- if(!insertVertex(v, worst, myFactory, allTets, vSizes)){
+ if(!insertVertex(v, worst, myFactory, allTets, vSizes,vSizesBGM)){
myFactory.changeTetRadius(allTets.begin(), 0.);
delete v;
}
diff --git a/Mesh/meshGRegionDelaunayInsertion.h b/Mesh/meshGRegionDelaunayInsertion.h
index 3377e514e59bec584e8465dc23dba6c9e2504593..104f4b0ae6b932a497f22ac0810860d50842a274 100644
--- a/Mesh/meshGRegionDelaunayInsertion.h
+++ b/Mesh/meshGRegionDelaunayInsertion.h
@@ -25,6 +25,7 @@
#include <map>
#include <stack>
#include "MElement.h"
+#include "BackgroundMesh.h"
#include "qualityMeasures.h"
//#define _GMSH_PRE_ALLOCATE_STRATEGY_ 1
@@ -85,7 +86,8 @@ class MTet4
double vol;
circum_radius = qmTet(t, qm, &vol);
}
- void setup(MTetrahedron *t, std::vector<double> &sizes)
+
+ void setup(MTetrahedron *t, std::vector<double> &sizes, std::vector<double> &sizesBGM)
{
base = t;
neigh[0] = neigh[1] = neigh[2] = neigh[3] = 0;
@@ -94,11 +96,18 @@ class MTet4
const double dx = base->getVertex(0)->x() - center[0];
const double dy = base->getVertex(0)->y() - center[1];
const double dz = base->getVertex(0)->z() - center[2];
- circum_radius = sqrt(dx * dx + dy * dy + dz * dz);
- double lc = 0.25 * (sizes[base->getVertex(0)->getNum()] +
- sizes[base->getVertex(1)->getNum()] +
- sizes[base->getVertex(2)->getNum()] +
- sizes[base->getVertex(3)->getNum()]);
+ circum_radius = sqrt ( dx*dx + dy*dy + dz*dz);
+
+ double lc1 = 0.25*(sizes [base->getVertex(0)->getNum()]+
+ sizes [base->getVertex(1)->getNum()]+
+ sizes [base->getVertex(2)->getNum()]+
+ sizes [base->getVertex(3)->getNum()]);
+ double lcBGM = 0.25*(sizesBGM [base->getVertex(0)->getNum()]+
+ sizesBGM [base->getVertex(1)->getNum()]+
+ sizesBGM [base->getVertex(2)->getNum()]+
+ sizesBGM [base->getVertex(3)->getNum()]);
+ double lc = Extend2dMeshIn3dVolumes() ? std::min(lc1, lcBGM) : lcBGM;
+
circum_radius /= lc;
deleted = false;
}
@@ -200,14 +209,14 @@ class MTet4Factory
delete [] allSlots;
#endif
}
- MTet4 *Create(MTetrahedron *t, std::vector<double> &sizes)
+ MTet4 * Create (MTetrahedron * t, std::vector<double> & sizes, std::vector<double> & sizesBGM)
{
#ifdef _GMSH_PRE_ALLOCATE_STRATEGY_
MTet4 *t4 = getAnEmptySlot();
#else
MTet4 *t4 = new MTet4;
#endif
- t4->setup(t, sizes);
+ t4->setup(t,sizes,sizesBGM);
return t4;
}
void Free(MTet4 *t)
diff --git a/Numeric/GaussLegendre1D.h b/Numeric/GaussLegendre1D.h
index 75513954642f4e29cfe60933a3f87d67e5a02389..a209e8a8fd579f4cca4665f80ed7e39272811e36 100644
--- a/Numeric/GaussLegendre1D.h
+++ b/Numeric/GaussLegendre1D.h
@@ -67,6 +67,84 @@ static double _GL_pt6[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};
+ /* 7 point rule points */
+ static double _GL_pt7[7]={
+ -9.491079123427585e-01,-7.415311855993945e-01,-4.058451513773972e-01, 0.000000000000000e+00, 4.058451513773972e-01, 7.415311855993945e-01, 9.491079123427585e-01};
+
+ /* 7 point rule weights */
+ static double _GL_wt7[7]={
+ 1.294849661688697e-01, 2.797053914892767e-01, 3.818300505051190e-01, 4.179591836734694e-01, 3.818300505051190e-01, 2.797053914892767e-01, 1.294849661688697e-01};
+
+ /* 8 point rule points */
+ static double _GL_pt8[8]={
+ -9.602898564975363e-01,-7.966664774136268e-01,-5.255324099163290e-01,-1.834346424956498e-01, 1.834346424956498e-01, 5.255324099163290e-01, 7.966664774136268e-01, 9.602898564975363e-01};
+
+ /* 8 point rule weights */
+ static double _GL_wt8[8]={
+ 1.012285362903768e-01, 2.223810344533745e-01, 3.137066458778874e-01, 3.626837833783620e-01, 3.626837833783620e-01, 3.137066458778874e-01, 2.223810344533745e-01, 1.012285362903768e-01};
+
+ /* 9 point rule points */
+ static double _GL_pt9[9]={
+ -9.681602395076261e-01,-8.360311073266359e-01,-6.133714327005905e-01,-3.242534234038089e-01, 0.000000000000000e+00, 3.242534234038089e-01, 6.133714327005905e-01, 8.360311073266359e-01, 9.681602395076261e-01};
+
+ /* 9 point rule weights */
+ static double _GL_wt9[9]={
+ 8.127438836157463e-02, 1.806481606948576e-01, 2.606106964029355e-01, 3.123470770400029e-01, 3.302393550012598e-01, 3.123470770400029e-01, 2.606106964029355e-01, 1.806481606948576e-01, 8.127438836157463e-02};
+
+ /* 10 point rule points */
+ static double _GL_pt10[10]={
+ -9.739065285171716e-01,-8.650633666889845e-01,-6.794095682990244e-01,-4.333953941292472e-01,-1.488743389816312e-01, 1.488743389816312e-01, 4.333953941292472e-01, 6.794095682990244e-01, 8.650633666889845e-01, 9.739065285171716e-01};
+
+ /* 10 point rule weights */
+ static double _GL_wt10[10]={
+ 6.667134430868774e-02, 1.494513491505805e-01, 2.190863625159822e-01, 2.692667193099962e-01, 2.955242247147529e-01, 2.955242247147529e-01, 2.692667193099962e-01, 2.190863625159822e-01, 1.494513491505805e-01, 6.667134430868774e-02};
+
+ /* 11 point rule points */
+ static double _GL_pt11[11]={
+ -9.782286581460570e-01,-8.870625997680953e-01,-7.301520055740494e-01,-5.190961292068118e-01,-2.695431559523450e-01, 0.000000000000000e+00, 2.695431559523450e-01, 5.190961292068118e-01, 7.301520055740494e-01, 8.870625997680953e-01, 9.782286581460570e-01};
+
+ /* 11 point rule weights */
+ static double _GL_wt11[11]={
+ 5.566856711617354e-02, 1.255803694649047e-01, 1.862902109277343e-01, 2.331937645919903e-01, 2.628045445102466e-01, 2.729250867779006e-01, 2.628045445102466e-01, 2.331937645919903e-01, 1.862902109277343e-01, 1.255803694649047e-01, 5.566856711617354e-02};
+
+ /* 12 point rule points */
+ static double _GL_pt12[12]={
+ -9.815606342467192e-01,-9.041172563704748e-01,-7.699026741943047e-01,-5.873179542866175e-01,-3.678314989981802e-01,-1.252334085114689e-01, 1.252334085114689e-01, 3.678314989981802e-01, 5.873179542866175e-01, 7.699026741943047e-01, 9.041172563704748e-01, 9.815606342467192e-01};
+
+ /* 12 point rule weights */
+ static double _GL_wt12[12]={
+ 4.717533638651183e-02, 1.069393259953182e-01, 1.600783285433463e-01, 2.031674267230658e-01, 2.334925365383548e-01, 2.491470458134029e-01, 2.491470458134029e-01, 2.334925365383548e-01, 2.031674267230658e-01, 1.600783285433463e-01, 1.069393259953182e-01, 4.717533638651183e-02};
+
+ /* 13 point rule points */
+ static double _GL_pt13[13]={
+ -9.841830547185881e-01,-9.175983992229780e-01,-8.015780907333099e-01,-6.423493394403402e-01,-4.484927510364468e-01,-2.304583159551348e-01, 1.232595164407831e-32, 2.304583159551348e-01, 4.484927510364468e-01, 6.423493394403402e-01, 8.015780907333099e-01, 9.175983992229780e-01, 9.841830547185881e-01};
+
+ /* 13 point rule weights */
+ static double _GL_wt13[13]={
+ 4.048400476531581e-02, 9.212149983772838e-02, 1.388735102197872e-01, 1.781459807619457e-01, 2.078160475368884e-01, 2.262831802628971e-01, 2.325515532308739e-01, 2.262831802628971e-01, 2.078160475368884e-01, 1.781459807619457e-01, 1.388735102197872e-01, 9.212149983772838e-02, 4.048400476531581e-02};
+
+ /* 14 point rule points */
+ static double _GL_pt14[14]={
+ -9.862838086968123e-01,-9.284348836635736e-01,-8.272013150697650e-01,-6.872929048116855e-01,-5.152486363581541e-01,-3.191123689278897e-01,-1.080549487073437e-01, 1.080549487073437e-01, 3.191123689278897e-01, 5.152486363581541e-01, 6.872929048116855e-01, 8.272013150697650e-01, 9.284348836635736e-01, 9.862838086968123e-01};
+
+ /* 14 point rule weights */
+ static double _GL_wt14[14]={
+ 3.511946033175199e-02, 8.015808715976037e-02, 1.215185706879031e-01, 1.572031671581936e-01, 1.855383974779378e-01, 2.051984637212955e-01, 2.152638534631578e-01, 2.152638534631578e-01, 2.051984637212955e-01, 1.855383974779378e-01, 1.572031671581936e-01, 1.215185706879031e-01, 8.015808715976037e-02, 3.511946033175199e-02};
+
+ /* 15 point rule points */
+ static double _GL_pt15[15]={
+ -9.879925180204854e-01,-9.372733924007060e-01,-8.482065834104272e-01,-7.244177313601701e-01,-5.709721726085388e-01,-3.941513470775634e-01,-2.011940939974345e-01, 1.232595164407831e-32, 2.011940939974345e-01, 3.941513470775634e-01, 5.709721726085388e-01, 7.244177313601701e-01, 8.482065834104272e-01, 9.372733924007060e-01, 9.879925180204854e-01};
+
+ /* 15 point rule weights */
+ static double _GL_wt15[15]={
+ 3.075324199611663e-02, 7.036604748810814e-02, 1.071592204671720e-01, 1.395706779261543e-01, 1.662692058169940e-01, 1.861610000155622e-01, 1.984314853271116e-01, 2.025782419255613e-01, 1.984314853271116e-01, 1.861610000155622e-01, 1.662692058169940e-01, 1.395706779261543e-01, 1.071592204671720e-01, 7.036604748810814e-02, 3.075324199611663e-02};
+/* 16 point rule points */
+ static double _GL_pt16[16]={
+ -9.894009349916499e-01,-9.445750230732326e-01,-8.656312023878318e-01,-7.554044083550030e-01,-6.178762444026438e-01,-4.580167776572274e-01,-2.816035507792589e-01,-9.501250983763744e-02, 9.501250983763744e-02, 2.816035507792589e-01, 4.580167776572274e-01, 6.178762444026438e-01, 7.554044083550030e-01, 8.656312023878318e-01, 9.445750230732326e-01, 9.894009349916499e-01};
+
+ /* 16 point rule weights */
+ static double _GL_wt16[16]={
+ 2.715245941175406e-02, 6.225352393864778e-02, 9.515851168249290e-02, 1.246289712555339e-01, 1.495959888165768e-01, 1.691565193950026e-01, 1.826034150449236e-01, 1.894506104550685e-01, 1.894506104550685e-01, 1.826034150449236e-01, 1.691565193950026e-01, 1.495959888165768e-01, 1.246289712555339e-01, 9.515851168249290e-02, 6.225352393864778e-02, 2.715245941175406e-02};
inline void gmshGaussLegendre1D(int nbQuadPoints, double **t, double **w)
{
@@ -95,6 +173,46 @@ inline void gmshGaussLegendre1D(int nbQuadPoints, double **t, double **w)
*t = _GL_pt6;
*w = _GL_wt6;
break;
+ case 7:
+ *t = _GL_pt7;
+ *w = _GL_wt7;
+ break;
+ case 8:
+ *t = _GL_pt8;
+ *w = _GL_wt8;
+ break;
+ case 9:
+ *t = _GL_pt9;
+ *w = _GL_wt9;
+ break;
+ case 10:
+ *t = _GL_pt10;
+ *w = _GL_wt10;
+ break;
+ case 11:
+ *t = _GL_pt11;
+ *w = _GL_wt11;
+ break;
+ case 12:
+ *t = _GL_pt12;
+ *w = _GL_wt12;
+ break;
+ case 13:
+ *t = _GL_pt13;
+ *w = _GL_wt13;
+ break;
+ case 14:
+ *t = _GL_pt14;
+ *w = _GL_wt14;
+ break;
+ case 15:
+ *t = _GL_pt15;
+ *w = _GL_wt15;
+ break;
+ case 16:
+ *t = _GL_pt16;
+ *w = _GL_wt16;
+ break;
default :
throw;
}