diff --git a/Mesh/Levy3D.cpp b/Mesh/Levy3D.cpp
index 0c81b133564736e5bee17bb4ca6a73ce8a525b2f..af47895f74227db64a2a08d3798553439c70e323 100755
--- a/Mesh/Levy3D.cpp
+++ b/Mesh/Levy3D.cpp
@@ -67,7 +67,7 @@ class LpCVT{
private:
std::vector<VoronoiElement> clipped;
fullMatrix<double> gauss_points;
- fullMatrix<double> gauss_weights;
+ fullVector<double> gauss_weights;
int gauss_num;
public:
LpCVT();
@@ -987,7 +987,7 @@ double LpCVT::F(VoronoiElement element,int p){
y = element.T(u,v,w,generator.y(),C1.y(),C2.y(),C3.y());
z = element.T(u,v,w,generator.z(),C1.z(),C2.z(),C3.z());
point = SPoint3(x,y,z);
- weight = gauss_weights(i,0);
+ weight = gauss_weights(i);
rho = h_to_rho(element.get_h(u,v,w),p); //h_to_rho(get_size(x,y,z),p);
energy = energy + weight*rho*f(generator,point,t,p);
}
@@ -1029,7 +1029,7 @@ SVector3 LpCVT::simple(VoronoiElement element,int p){
y = element.T(u,v,w,generator.y(),C1.y(),C2.y(),C3.y());
z = element.T(u,v,w,generator.z(),C1.z(),C2.z(),C3.z());
point = SPoint3(x,y,z);
- weight = gauss_weights(i,0);
+ weight = gauss_weights(i);
rho = h_to_rho(element.get_h(u,v,w),p); //h_to_rho(get_size(x,y,z),p);
comp_x = comp_x + weight*rho*df_dx(generator,point,t,p);
comp_y = comp_y + weight*rho*df_dy(generator,point,t,p);
@@ -1081,7 +1081,7 @@ SVector3 LpCVT::dF_dC1(VoronoiElement element,int p){
y = element.T(u,v,w,gy,C1.y(),C2.y(),C3.y());
z = element.T(u,v,w,gz,C1.z(),C2.z(),C3.z());
point = SPoint3(x,y,z);
- weight = gauss_weights(i,0);
+ weight = gauss_weights(i);
rho = h_to_rho(element.get_h(u,v,w),p); //h_to_rho(get_size(x,y,z),p);
drho_dx = (-p-3)*pow_int(1.0/element.get_h(u,v,w),p+4)*element.get_dh_dx(); //get_drho_dx(x,y,z,p);
drho_dy = (-p-3)*pow_int(1.0/element.get_h(u,v,w),p+4)*element.get_dh_dy(); //get_drho_dy(x,y,z,p);
@@ -1140,7 +1140,7 @@ SVector3 LpCVT::dF_dC2(VoronoiElement element,int p){
y = element.T(u,v,w,generator.y(),C1.y(),C2.y(),C3.y());
z = element.T(u,v,w,generator.z(),C1.z(),C2.z(),C3.z());
point = SPoint3(x,y,z);
- weight = gauss_weights(i,0);
+ weight = gauss_weights(i);
rho = h_to_rho(element.get_h(u,v,w),p); //h_to_rho(get_size(x,y,z),p);
drho_dx = (-p-3)*pow_int(1.0/element.get_h(u,v,w),p+4)*element.get_dh_dx(); //get_drho_dx(x,y,z,p);
drho_dy = (-p-3)*pow_int(1.0/element.get_h(u,v,w),p+4)*element.get_dh_dy(); //get_drho_dy(x,y,z,p);
@@ -1199,7 +1199,7 @@ SVector3 LpCVT::dF_dC3(VoronoiElement element,int p){
y = element.T(u,v,w,gy,C1.y(),C2.y(),C3.y());
z = element.T(u,v,w,gz,C1.z(),C2.z(),C3.z());
point = SPoint3(x,y,z);
- weight = gauss_weights(i,0);
+ weight = gauss_weights(i);
rho = h_to_rho(element.get_h(u,v,w),p); //h_to_rho(get_size(x,y,z),p);
drho_dx = (-p-3)*pow_int(1.0/element.get_h(u,v,w),p+4)*element.get_dh_dx(); //get_drho_dx(x,y,z,p);
drho_dy = (-p-3)*pow_int(1.0/element.get_h(u,v,w),p+4)*element.get_dh_dy(); //get_drho_dy(x,y,z,p);
diff --git a/Mesh/meshGFaceLloyd.cpp b/Mesh/meshGFaceLloyd.cpp
index 8c70482d3376d60ec8828189fb9001c87a8ceca4..5c3c288b2d4ab6f1c90ef1299a4a35f2580a5730 100644
--- a/Mesh/meshGFaceLloyd.cpp
+++ b/Mesh/meshGFaceLloyd.cpp
@@ -189,7 +189,7 @@ class lpcvt{
std::vector<double> angles;
std::vector<voronoi_cell> temp;
fullMatrix<double> gauss_points;
- fullMatrix<double> gauss_weights;
+ fullVector<double> gauss_weights;
int gauss_num;
public :
lpcvt();
@@ -1440,7 +1440,7 @@ double lpcvt::F(voronoi_element element,int p){
x = Tx(u,v,generator,C1,C2);
y = Ty(u,v,generator,C1,C2);
point = SPoint2(x,y);
- weight = gauss_weights(i,0);
+ weight = gauss_weights(i);
rho = element.get_rho(u,v);
energy = energy + weight*rho*f(generator,point,m,p);
}
@@ -1480,7 +1480,7 @@ SVector3 lpcvt::simple(voronoi_element element,int p){
x = Tx(u,v,generator,C1,C2);
y = Ty(u,v,generator,C1,C2);
point = SPoint2(x,y);
- weight = gauss_weights(i,0);
+ weight = gauss_weights(i);
rho = element.get_rho(u,v);
comp_x = comp_x + weight*rho*df_dx(generator,point,m,p);
comp_y = comp_y + weight*rho*df_dy(generator,point,m,p);
@@ -1525,7 +1525,7 @@ SVector3 lpcvt::dF_dC1(voronoi_element element,int p){
x = Tx(u,v,generator,C1,C2);
y = Ty(u,v,generator,C1,C2);
point = SPoint2(x,y);
- weight = gauss_weights(i,0);
+ weight = gauss_weights(i);
rho = element.get_rho(u,v);
drho_dx = element.get_drho_dx();
drho_dy = element.get_drho_dy();
@@ -1575,7 +1575,7 @@ SVector3 lpcvt::dF_dC2(voronoi_element element,int p){
x = Tx(u,v,generator,C1,C2);
y = Ty(u,v,generator,C1,C2);
point = SPoint2(x,y);
- weight = gauss_weights(i,0);
+ weight = gauss_weights(i);
rho = element.get_rho(u,v);
drho_dx = element.get_drho_dx();
drho_dy = element.get_drho_dy();
diff --git a/Numeric/GaussIntegration.cpp b/Numeric/GaussIntegration.cpp
index 61531226ced6ca1f952f8f34a254d32d833bad63..756ca222d2f73a9c3796cd4674feffbebd1a6034 100644
--- a/Numeric/GaussIntegration.cpp
+++ b/Numeric/GaussIntegration.cpp
@@ -7,7 +7,7 @@
#include "GmshDefines.h"
static void pts2fullMatrix(int npts, IntPt *pts, fullMatrix<double> &pMatrix,
- fullMatrix<double> &wMatrix)
+ fullVector<double> &wMatrix)
{
pMatrix.resize(npts,3);
wMatrix.resize(npts,1);
@@ -15,48 +15,48 @@ static void pts2fullMatrix(int npts, IntPt *pts, fullMatrix<double> &pMatrix,
pMatrix(i, 0) = pts[i].pt[0];
pMatrix(i, 1) = pts[i].pt[1];
pMatrix(i, 2) = pts[i].pt[2];
- wMatrix(i, 0) = pts[i].weight;
+ wMatrix(i) = pts[i].weight;
}
}
void gaussIntegration::getTriangle(int order, fullMatrix<double> &pts,
- fullMatrix<double> &weights)
+ fullVector<double> &weights)
{
pts2fullMatrix(getNGQTPts(order),getGQTPts(order),pts,weights);
}
void gaussIntegration::getLine(int order, fullMatrix<double> &pts,
- fullMatrix<double> &weights)
+ fullVector<double> &weights)
{
pts2fullMatrix(getNGQLPts(order),getGQLPts(order),pts,weights);
}
void gaussIntegration::getQuad(int order, fullMatrix<double> &pts,
- fullMatrix<double> &weights)
+ fullVector<double> &weights)
{
pts2fullMatrix(getNGQQPts(order),getGQQPts(order),pts,weights);
}
void gaussIntegration::getTetrahedron(int order, fullMatrix<double> &pts,
- fullMatrix<double> &weights)
+ fullVector<double> &weights)
{
pts2fullMatrix(getNGQTetPts(order),getGQTetPts(order),pts,weights);
}
void gaussIntegration::getHexahedron(int order, fullMatrix<double> &pts,
- fullMatrix<double> &weights)
+ fullVector<double> &weights)
{
pts2fullMatrix(getNGQHPts(order),getGQHPts(order),pts,weights);
}
void gaussIntegration::getPrism(int order, fullMatrix<double> &pts,
- fullMatrix<double> &weights)
+ fullVector<double> &weights)
{
pts2fullMatrix(getNGQPriPts(order),getGQPriPts(order),pts,weights);
}
void gaussIntegration::get(int elementType, int order, fullMatrix<double> &pts,
- fullMatrix<double> &weights)
+ fullVector<double> &weights)
{
switch (elementType) {
case TYPE_TRI : pts2fullMatrix(getNGQTPts(order), getGQTPts(order), pts, weights); break;
@@ -67,7 +67,7 @@ void gaussIntegration::get(int elementType, int order, fullMatrix<double> &pts,
case TYPE_PRI : pts2fullMatrix(getNGQPriPts(order), getGQPriPts(order), pts, weights); break;
case TYPE_PNT :
weights.resize(1,1);
- weights(0,0) = 1.;
+ weights(0) = 1.;
pts.resize(1,3);
break;
default : Msg::Error("No integration rules defined for type %i", elementType);
diff --git a/Numeric/GaussIntegration.h b/Numeric/GaussIntegration.h
index a3d5ca427ca385af2ed857ae8119131088b522b4..9c9079965a42016e81752c577968e3a2a4df9f2a 100644
--- a/Numeric/GaussIntegration.h
+++ b/Numeric/GaussIntegration.h
@@ -37,19 +37,19 @@ IntPt *getGQHPts(int order);
class gaussIntegration {
public:
static void get(int elementType, int order, fullMatrix<double> &pts,
- fullMatrix<double> &weights);
+ fullVector<double> &weights);
static void getTriangle(int order, fullMatrix<double> &pts,
- fullMatrix<double> &weights);
+ fullVector<double> &weights);
static void getLine(int order, fullMatrix<double> &pts,
- fullMatrix<double> &weights);
+ fullVector<double> &weights);
static void getQuad(int order, fullMatrix<double> &pts,
- fullMatrix<double> &weights);
+ fullVector<double> &weights);
static void getTetrahedron(int order, fullMatrix<double> &pts,
- fullMatrix<double> &weights);
+ fullVector<double> &weights);
static void getHexahedron(int order, fullMatrix<double> &pts,
- fullMatrix<double> &weights);
+ fullVector<double> &weights);
static void getPrism(int order, fullMatrix<double> &pts,
- fullMatrix<double> &weights);
+ fullVector<double> &weights);
};
#endif
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index af521711d35309f95d58a2d0e7d90bb56b4907b4..9c854dd6ac88a8db686a4b85b9c2b86252092171 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -1247,6 +1247,7 @@ static void getFaceClosurePrism(int iFace, int iSign, int iRotate,
bool isTriangle = iFace<2;
int nNodes = isTriangle ? (order+1)*(order+2)/2 : (order+1)*(order+1);
closure.clear();
+ if (isTriangle && iRotate > 2) return;
closure.resize(nNodes);
if (order==0) {
closure[0] = 0;
@@ -1301,6 +1302,8 @@ static void generateFaceClosurePrismFull(polynomialBasis::clCont &closureFull,
for (unsigned int i = 0; i < closure.size(); i++) {
std::vector<int> &clFull = closureFull[i];
std::vector<int> &cl = closure[i];
+ if (cl.size() == 0)
+ continue;
clFull.resize(6, -1);
int &ref = cl.size() == 3 ? ref3 : (cl[0] / 3 + cl[1] / 3) % 2 ? ref4b : ref4a;
if (ref == -1)
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index 64071e0b97328346d2d8db0904e270562f554eaf..f0f42c2ddb4ab480af087ca1cbdaf827329027dd 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -96,6 +96,10 @@ class polynomialBasis
inline int getNumShapeFunctions() const {return coefficients.size1();}
// for a given face/edge, with both a sign and a rotation, give an
// ordered list of nodes on this face/edge
+ inline int getClosureType(int id) const
+ {
+ return closures[id].type;
+ }
inline const std::vector<int> &getClosure(int id) const
{
return closures[id];