From b827fb021cf787d0d756eab24ce9a49743fdfeb4 Mon Sep 17 00:00:00 2001
From: Amaury Johnan <amjohnen@gmail.com>
Date: Mon, 10 Jun 2013 14:07:00 +0000
Subject: [PATCH] Create gmshGenerateMonomials. Goal: at the end, have only one
algo for generation of Points/Monomials/Exponents/.. per element
---
Geo/MQuadrangle.h | 2 +-
Geo/MTetrahedron.h | 34 +-
Numeric/nodalBasis.cpp | 3 +-
Numeric/nodalBasis.h | 39 +--
Numeric/pointsGenerators.cpp | 310 ++++++++++++++++++
Numeric/pointsGenerators.h | 32 +-
Numeric/polynomialBasis.h | 43 +--
.../voro++/examples/degenerate/degenerate2.cc | 4 +-
8 files changed, 369 insertions(+), 98 deletions(-)
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index 19591a4342..14e3bf691e 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -343,7 +343,7 @@ class MQuadrangle9 : public MQuadrangle {
* | | N = total number of vertices
* 4+3E 3+2E
* | | Interior vertex numbers
- * ... 4+4e to N-1 ... for edge 0 <= i <= 3: 4+i*E to 3+(i+1)*E
+ * ... 4+4E to N-1 ... for edge 0 <= i <= 3: 4+i*E to 3+(i+1)*E
* | | in volume : 4+4*E to N-1
* 3+4E 4+E
* | |
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index 5c6a7b15db..920a88bc2f 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -275,37 +275,19 @@ class MTetrahedron10 : public MTetrahedron {
}
};
-/*
- * MTetrahedronN FIXME: check the plot
- *
- * 2
- * ,/|`\
- * ,/ | `\ E = order - 1
- * ,/ '. `\ C = 4 + 6*E
- * ,/ | `\ F = ((order - 1)*(order - 2))/2
- * ,/ | `\ N = total number of vertices
- * 0-----------'.--------1
- * `\. | ,/ Interior vertex numbers
- * `\. | ,/ for edge 0 <= i <= 5: 4+i*E to 3+(i+1)*E
- * `\. '. ,/ for face 0 <= j <= 3: C+j*F to C-1+(j+1)*F
- * `\. |/ in volume : C+4*F to N-1
- * `3
+/* tet order 3 FIXME: check the plot
*
- */
-
-/* tet order 3
-
* 2
* ,/|`\
- * ,5 | `6 E = order - 1
- * ,/ 12 `\ C = 4 + 6*E
- * ,4 | `7 F = ((order - 1)*(order - 2))/2
+ * ,8 | `7 E = order - 1
+ * ,/ 13 `\ C = 4 + 6*E
+ * ,9 | `6 F = ((order - 1)*(order - 2))/2
* ,/ | `\ N = total number of vertices
- * 0-----9-----'.--8-----1
+ * 0-----4-----'.--5-----1
* `\. | ,/ Interior vertex numbers
- * 10. 13 ,14 for edge 0 <= i <= 5: 4+i*E to 3+(i+1)*E
- * `\. '. 15 for face 0 <= j <= 3: C+j*F to C-1+(j+1)*F
- * 11\.|/ in volume : C+4*F to N-1
+ * 11. 12 ,15 for edge 0 <= i <= 5: 4+i*E to 4+(i+1)*E-1
+ * `\. '. 14 for face 0 <= j <= 3: C+j*F to C+(j+1)*F-1
+ * 10\.|/ in volume : C+4*F to N-1
* `3
*
*/
diff --git a/Numeric/nodalBasis.cpp b/Numeric/nodalBasis.cpp
index 38126f89c9..8a2c1b3c65 100644
--- a/Numeric/nodalBasis.cpp
+++ b/Numeric/nodalBasis.cpp
@@ -1372,9 +1372,8 @@ static void generateClosureOrder0(nodalBasis::clCont &closure, int nb)
nodalBasis::nodalBasis(int tag)
{
-
type = tag;
-
+
switch (tag) {
case MSH_PNT : parentType = TYPE_PNT; order = 0; serendip = false; break;
case MSH_LIN_1 : parentType = TYPE_LIN; order = 0; serendip = false; break;
diff --git a/Numeric/nodalBasis.h b/Numeric/nodalBasis.h
index 5bc687fdcb..eb545bbeaa 100644
--- a/Numeric/nodalBasis.h
+++ b/Numeric/nodalBasis.h
@@ -10,33 +10,26 @@
#include "GmshDefines.h"
class nodalBasis {
-
public:
int type, parentType, order, dimension, numFaces;
bool serendip;
fullMatrix<double> points;
-
+
nodalBasis(int tag);
virtual ~nodalBasis() {}
-
- // Basis functions evaluation
- virtual void f(double u, double v, double w, double *sf) const {Msg::Fatal("Not implemented");};
- virtual void f(const fullMatrix<double> &coord, fullMatrix<double> &sf) const {Msg::Fatal("Not implemented");};
-
- // Basis functions gradients evaluation
- virtual void df(double u, double v, double w, double grads[][3]) const {Msg::Fatal("Not implemented");};
- virtual void df(const fullMatrix<double> &coord, fullMatrix<double> &dfm) const {Msg::Fatal("Not implemented");};
+ virtual int getNumShapeFunctions() const = 0;
+
+ // Basis functions & gradients evaluation
+ virtual void f(double u, double v, double w, double *sf) const = 0;
+ virtual void f(const fullMatrix<double> &coord, fullMatrix<double> &sf) const = 0;
+ virtual void df(double u, double v, double w, double grads[][3]) const = 0;
+ virtual void df(const fullMatrix<double> &coord, fullMatrix<double> &dfm) const = 0;
virtual void ddf(double u, double v, double w, double grads[][3][3]) const {Msg::Fatal("Not implemented");};
virtual void dddf(double u, double v, double w, double grads[][3][3][3]) const {Msg::Fatal("Not implemented");};
-
- virtual int getNumShapeFunctions() const {Msg::Fatal("Not implemented"); return -1;}
-
- class closure : public std::vector<int> {
- public:
- int type;
- };
- typedef std::vector<closure> clCont;
+
+
+
// closures is the list of the nodes of each face, in the order of
// the polynomialBasis of the face; fullClosures is mapping of the
// nodes of the element that rotates the element so that the
@@ -45,9 +38,14 @@ class nodalBasis {
// fullCLosure[i] rotates the element so that the considered face
// becomes the closureRef[i]-th face (the first tringle or the first
// quad face)
+ class closure : public std::vector<int> {
+ public:
+ int type;
+ };
+ typedef std::vector<closure> clCont;
clCont closures, fullClosures;
std::vector<int> closureRef;
-
+
// for a given face/edge, with both a sign and a rotation, give an
// ordered list of nodes on this face/edge
virtual int getClosureType(int id) const { return closures[id].type; }
@@ -55,9 +53,8 @@ class nodalBasis {
virtual const std::vector<int> &getFullClosure(int id) const { return fullClosures[id]; }
inline int getClosureId(int iFace, int iSign=1, int iRot=0) const;
inline void breakClosureId(int i, int &iFace, int &iSign, int &iRot) const;
-
+
static inline int getTag(int parentTag, int order, bool serendip = false);
-
};
diff --git a/Numeric/pointsGenerators.cpp b/Numeric/pointsGenerators.cpp
index 6c3db2a1a7..5b4e4221ac 100644
--- a/Numeric/pointsGenerators.cpp
+++ b/Numeric/pointsGenerators.cpp
@@ -4,6 +4,10 @@
// bugs and problems to the public mailing list <gmsh@geuz.org>.
#include "pointsGenerators.h"
+#include "MTriangle.h"
+#include "MQuadrangle.h"
+#include "MTetrahedron.h"
+#include "MHexahedron.h"
/* --- Lines --- */
@@ -779,3 +783,309 @@ fullMatrix<double> gmshGeneratePointsPyramid(int order, bool serendip)
}
return points;
}
+
+
+
+
+fullMatrix<int> gmshGenerateMonomialsLine(int order)
+{
+ fullMatrix<int> monomials(order + 1, 1);
+ monomials(0,0) = 0;
+ if (order > 0) {
+ monomials(1, 0) = order;
+ for (int i = 2; i < order + 1; i++) monomials(i, 0) = i-1;
+ }
+ return monomials;
+}
+fullMatrix<int> gmshGenerateMonomialsTriangle(int order, bool serendip)
+{
+ int nbMonomials = serendip ? 3 * order : (order + 1) * (order + 2) / 2;
+ fullMatrix<int> monomials(nbMonomials, 2);
+
+ monomials(0, 0) = 0;
+ monomials(0, 1) = 0;
+
+ if (order > 0) {
+ monomials(1, 0) = order;
+ monomials(1, 1) = 0;
+
+ monomials(2, 0) = 0;
+ monomials(2, 1) = order;
+
+ if (order > 1) {
+ int index = 3;
+ for (int iedge = 0; iedge < 3; ++iedge) {
+ int i0 = MTriangle::edges_tri(iedge, 0);
+ int i1 = MTriangle::edges_tri(iedge, 1);
+
+ int u_0 = (monomials(i1,0)-monomials(i0,0)) / order;
+ int u_1 = (monomials(i1,1)-monomials(i0,1)) / order;
+
+ for (int i = 1; i < order; ++i, ++index) {
+ monomials(index, 0) = monomials(i0, 0) + u_0 * i;
+ monomials(index, 1) = monomials(i0, 1) + u_1 * i;
+ }
+ }
+ fullMatrix<int> inner = gmshGenerateMonomialsQuadrangle(order-2);
+ inner.add(1);
+ monomials.copy(inner, 0, nbMonomials - index, 0, 2, index, 0);
+ }
+ }
+ return monomials;
+}
+fullMatrix<int> gmshGenerateMonomialsQuadrangle(int order)
+{
+ int nbMonomials = (order+1)*(order+1);
+ fullMatrix<int> monomials(nbMonomials, 2);
+
+ monomials(0, 0) = 0;
+ monomials(0, 1) = 0;
+
+ if (order > 0) {
+ monomials(1, 0) = order;
+ monomials(1, 1) = 0;
+
+ monomials(2, 0) = order;
+ monomials(2, 1) = order;
+
+ monomials(3, 0) = 0;
+ monomials(3, 1) = order;
+
+ if (order > 1) {
+ int index = 4;
+ for (int iedge = 0; iedge < 4; ++iedge) {
+ int i0 = MQuadrangle::edges_quad(iedge, 0);
+ int i1 = MQuadrangle::edges_quad(iedge, 1);
+
+ int u_0 = (monomials(i1,0)-monomials(i0,0)) / order;
+ int u_1 = (monomials(i1,1)-monomials(i0,1)) / order;
+
+ for (int i = 1; i < order; ++i, ++index) {
+ monomials(index, 0) = monomials(i0, 0) + u_0 * i;
+ monomials(index, 1) = monomials(i0, 1) + u_1 * i;
+ }
+ }
+ fullMatrix<int> inner = gmshGenerateMonomialsQuadrangle(order-2);
+ inner.add(1);
+ monomials.copy(inner, 0, nbMonomials - index, 0, 2, index, 0);
+ }
+ }
+ return monomials;
+}
+//KH : caveat : node coordinates are not yet coherent with node numbering associated
+// to numbering of principal vertices of face !!!!
+fullMatrix<int> gmshGenerateMonomialsTetrahedron(int order, bool serendip)
+{
+ int nbMonomials =
+ (serendip ?
+ 4 + 6 * std::max(0, order - 1) + 4 * std::max(0, (order - 2) * (order - 1) / 2) :
+ (order + 1) * (order + 2) * (order + 3) / 6);
+
+ fullMatrix<int> monomials(nbMonomials, 3);
+
+ monomials(0, 0) = 0;
+ monomials(0, 1) = 0;
+ monomials(0, 2) = 0;
+
+ if (order > 0) {
+ monomials(1, 0) = order;
+ monomials(1, 1) = 0;
+ monomials(1, 2) = 0;
+
+ monomials(2, 0) = 0;
+ monomials(2, 1) = order;
+ monomials(2, 2) = 0;
+
+ monomials(3, 0) = 0;
+ monomials(3, 1) = 0;
+ monomials(3, 2) = order;
+
+ // the template has been defined in table edges_tetra and faces_tetra (MElement.h)
+
+ if (order > 1) {
+ int index = 4;
+ for (int iedge = 0; iedge < 6; ++iedge) {
+ int i0 = MTetrahedron::edges_tetra(iedge, 0);
+ int i1 = MTetrahedron::edges_tetra(iedge, 1);
+
+ int u[3];
+ u[0] = (monomials(i1,0)-monomials(i0,0)) / order;
+ u[1] = (monomials(i1,1)-monomials(i0,1)) / order;
+ u[2] = (monomials(i1,2)-monomials(i0,2)) / order;
+
+ for (int i = 1; i < order; ++i, ++index) {
+ monomials(index, 0) = monomials(i0, 0) + u[0] * i;
+ monomials(index, 1) = monomials(i0, 1) + u[1] * i;
+ monomials(index, 2) = monomials(i0, 2) + u[2] * i;
+ }
+ }
+
+ if (order > 2) {
+ fullMatrix<int> dudv = gmshGenerateMonomialsTriangle(order - 3);
+
+ for (int iface = 0; iface < 4; ++iface) {
+ int i0 = MTetrahedron::faces_tetra(iface, 0);
+ int i1 = MTetrahedron::faces_tetra(iface, 1);
+ int i2 = MTetrahedron::faces_tetra(iface, 2);
+
+ int u[3];
+ u[0] = (monomials(i1, 0) - monomials(i0, 0)) / order;
+ u[1] = (monomials(i1, 1) - monomials(i0, 1)) / order;
+ u[2] = (monomials(i1, 2) - monomials(i0, 2)) / order;
+ int v[3];
+ v[0] = (monomials(i2, 0) - monomials(i0, 0)) / order;
+ v[1] = (monomials(i2, 1) - monomials(i0, 1)) / order;
+ v[2] = (monomials(i2, 2) - monomials(i0, 2)) / order;
+
+ for (int i = 0; i < dudv.size1(); ++i, ++index) {
+ monomials(index, 0) = monomials(i0, 0) + u[0] * dudv(i, 0) + v[0] * dudv(i, 1);
+ monomials(index, 1) = monomials(i0, 1) + u[1] * dudv(i, 0) + v[1] * dudv(i, 1);
+ monomials(index, 2) = monomials(i0, 2) + u[2] * dudv(i, 0) + v[2] * dudv(i, 1);
+ }
+ }
+
+ if (!serendip && order > 3) {
+ fullMatrix<int> inner = gmshGenerateMonomialsTetrahedron(order - 4);
+ inner.add(1);
+ monomials.copy(inner, 0, nbMonomials - index, 0, 3, index, 0);
+ }
+ }
+ }
+ }
+ return monomials;
+}
+/*
+fullMatrix<int> gmshGenerateMonomialsPrism(int order)
+{
+ const int prism18Pts[18][3] = {
+ {0, 0, 0}, // 0
+ {2, 0, 0}, // 1
+ {0, 2, 0}, // 2
+ {0, 0, 2}, // 3
+ {2, 0, 2}, // 4
+ {0, 2, 2}, // 5
+ {1, 0, 0}, // 6
+ {0, 1, 0}, // 7
+ {0, 0, 1}, // 8
+ {1, 1, 0}, // 9
+ {2, 0, 1}, // 10
+ {0, 2, 1}, // 11
+ {1, 0, 2}, // 12
+ {0, 1, 2}, // 13
+ {1, 1, 2}, // 14
+ {1, 0, 1}, // 15
+ {0, 1, 1}, // 16
+ {1, 1, 1}, // 17
+ };
+
+ int nbMonomials = (order + 1)*(order + 1)*(order + 2)/2;
+ fullMatrix<int> monomials(nbMonomials, 3);
+
+ int index = 0;
+ fullMatrix<int> triMonomials = gmshGenerateMonomialsTriangle(order,false);
+ fullMatrix<int> lineMonomials = gmshGenerateMonomialsLine(order);
+
+ if (order == 2)
+ for (int i =0; i<18; i++)
+ for (int j=0; j<3;j++)
+ monomials(i,j) = prism18Pts[i][j];
+ else
+ for (int j = 0; j <lineMonomials.size1() ; j++) {
+ for (int i = 0; i < triMonomials.size1(); i++) {
+ monomials(index,0) = triMonomials(i,0);
+ monomials(index,1) = triMonomials(i,1);
+ monomials(index,2) = lineMonomials(j,0);
+ index++;
+ }
+ }
+
+ return monomials;
+
+}
+*/
+fullMatrix<int> gmshGenerateMonomialsHexahedron(int order)
+{
+ int nbMonomials = (order+1)*(order+1)*(order+1);
+ fullMatrix<int> monomials(nbMonomials, 3);
+
+ monomials(0, 0) = 0;
+ monomials(0, 1) = 0;
+ monomials(0, 2) = 0;
+
+ if (order > 0) {
+ monomials(1, 0) = order;
+ monomials(1, 1) = 0;
+ monomials(1, 2) = 0;
+
+ monomials(2, 0) = order;
+ monomials(2, 1) = order;
+ monomials(2, 2) = 0;
+
+ monomials(3, 0) = 0;
+ monomials(3, 1) = order;
+ monomials(3, 2) = 0;
+
+ monomials(4, 0) = 0;
+ monomials(4, 1) = 0;
+ monomials(4, 2) = order;
+
+ monomials(5, 0) = order;
+ monomials(5, 1) = 0;
+ monomials(5, 2) = order;
+
+ monomials(6, 0) = order;
+ monomials(6, 1) = order;
+ monomials(6, 2) = order;
+
+ monomials(7, 0) = 0;
+ monomials(7, 1) = order;
+ monomials(7, 2) = order;
+
+ if (order > 1) {
+ int index = 8;
+ for (int iedge = 0; iedge < 12; ++iedge) {
+ int i0 = MHexahedron::edges_hexa(iedge, 0);
+ int i1 = MHexahedron::edges_hexa(iedge, 1);
+
+ int u_1 = (monomials(i1,0)-monomials(i0,0)) / order;
+ int u_2 = (monomials(i1,1)-monomials(i0,1)) / order;
+ int u_3 = (monomials(i1,2)-monomials(i0,2)) / order;
+
+ for (int i = 1; i < order; ++i, ++index) {
+ monomials(index, 0) = monomials(i0, 0) + i * u_1;
+ monomials(index, 1) = monomials(i0, 1) + i * u_2;
+ monomials(index, 2) = monomials(i0, 2) + i * u_3;
+ }
+ }
+
+ fullMatrix<int> dudv = gmshGenerateMonomialsQuadrangle(order - 3);
+
+ for (int iface = 0; iface < 6; ++iface) {
+ int i0 = MHexahedron::faces_hexa(iface, 0);
+ int i1 = MHexahedron::faces_hexa(iface, 1);
+ int i3 = MHexahedron::faces_hexa(iface, 3);
+
+ int u[3];
+ u[0] = (monomials(i1, 0) - monomials(i0, 0)) / order;
+ u[1] = (monomials(i1, 1) - monomials(i0, 1)) / order;
+ u[2] = (monomials(i1, 2) - monomials(i0, 2)) / order;
+ int v[3];
+ v[0] = (monomials(i3, 0) - monomials(i0, 0)) / order;
+ v[1] = (monomials(i3, 1) - monomials(i0, 1)) / order;
+ v[2] = (monomials(i3, 2) - monomials(i0, 2)) / order;
+
+ for (int i = 0; i < dudv.size1(); ++i, ++index) {
+ monomials(index, 0) = monomials(i0, 0) + u[0] * dudv(i, 0) + v[0] * dudv(i, 1);
+ monomials(index, 1) = monomials(i0, 1) + u[1] * dudv(i, 0) + v[1] * dudv(i, 1);
+ monomials(index, 2) = monomials(i0, 2) + u[2] * dudv(i, 0) + v[2] * dudv(i, 1);
+ }
+ }
+
+ fullMatrix<int> inner = gmshGenerateMonomialsHexahedron(order - 2);
+ inner.add(1);
+ monomials.copy(inner, 0, nbMonomials - index, 0, 3, index, 0);
+ }
+ }
+ return monomials;
+}
+
diff --git a/Numeric/pointsGenerators.h b/Numeric/pointsGenerators.h
index 583b874d15..5aa57b2c7f 100644
--- a/Numeric/pointsGenerators.h
+++ b/Numeric/pointsGenerators.h
@@ -11,34 +11,36 @@
/*
* Functions to generate point distributions on
* the references elements, for all orders.
+ * &
+ * Functions generating exponents of Pascal monomials
+ * in the same order than Gmsh Points.
*/
-/* --- Lines --- */
-
+// Points
fullMatrix<double> gmshGeneratePointsLine(int order);
-/* --- Triangles --- */
-
fullMatrix<double> gmshGeneratePointsTriangle(int order, bool serendip);
-
-/* --- Quadrangles --- */
-
fullMatrix<double> gmshGeneratePointsQuadrangle(int order, bool serendip);
-/* --- Tetahedra --- */
-
fullMatrix<double> gmshGeneratePointsTetrahedron(int order, bool serendip);
+fullMatrix<double> gmshGeneratePointsHexahedron(int order, bool serendip);
+fullMatrix<double> gmshGeneratePointsPrism(int order, bool serendip);
-/* --- Hexahedra --- */
+fullMatrix<double> gmshGeneratePointsPyramid(int order, bool serendip);
-fullMatrix<double> gmshGeneratePointsHexahedron(int order, bool serendip);
+// Monomial exponents
+fullMatrix<int> gmshGenerateMonomialsLine(int order);
-/* --- Prisms --- */
+fullMatrix<int> gmshGenerateMonomialsTriangle(int order, bool serendip = false);
+fullMatrix<int> gmshGenerateMonomialsQuadrangle(int order);
+//fullMatrix<int> gmshGenerateMonomialsQuadSerendipity(int order);
-fullMatrix<double> gmshGeneratePointsPrism(int order, bool serendip);
+fullMatrix<int> gmshGenerateMonomialsTetrahedron(int order, bool serendip = false);
+fullMatrix<int> gmshGenerateMonomialsHexahedron(int order);
+//fullMatrix<int> gmshGenerateMonomialsPrism(int order);
+
+//fullMatrix<int> gmshGenerateMonomialsPyramid(int order, bool serendip);
-/* --- Pyramids --- */
-fullMatrix<double> gmshGeneratePointsPyramid(int order, bool serendip);
#endif
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index 0c09ef633c..395f55ada3 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -64,52 +64,33 @@ inline double pow_int(const double &a, const int &n)
}
}
-fullMatrix<double> generate1DMonomials(int order);
-
-
-
-
class polynomialBasis : public nodalBasis
{
- // integrationOrder, closureId => df/dXi
-// mutable std::map<int,std::vector<fullMatrix<double> > > _dfAtFace;
public:
// for now the only implemented polynomial basis are nodal poly
// basis, we use the type of the corresponding gmsh element as type
- //int type, parentType, order, dimension;
- //bool serendip;
-
- //fullMatrix<double> points;
fullMatrix<double> monomials;
fullMatrix<double> coefficients;
-
+
polynomialBasis(int tag);
~polynomialBasis();
-
+
+ virtual int getNumShapeFunctions() const;
+
virtual void f(double u, double v, double w, double *sf) const;
virtual void f(const fullMatrix<double> &coord, fullMatrix<double> &sf) const;
virtual void df(const fullMatrix<double> &coord, fullMatrix<double> &dfm) const;
virtual void df(double u, double v, double w, double grads[][3]) const;
virtual void ddf(double u, double v, double w, double hess[][3][3]) const;
virtual void dddf(double u, double v, double w, double third[][3][3][3]) const;
-
- virtual int getNumShapeFunctions() const;
-
- inline void evaluateMonomials(double u, double v, double w, double p[]) const;
-
-};
-
-
-
-inline void polynomialBasis::evaluateMonomials(double u, double v, double w, double p[]) const
-{
- for (int j = 0; j < monomials.size1(); j++) {
- p[j] = pow_int(u, (int)monomials(j, 0));
- if (monomials.size2() > 1) p[j] *= pow_int(v, (int)monomials(j, 1));
- if (monomials.size2() > 2) p[j] *= pow_int(w, (int)monomials(j, 2));
+
+ inline void evaluateMonomials(double u, double v, double w, double p[]) const {
+ for (int j = 0; j < monomials.size1(); j++) {
+ p[j] = pow_int(u, (int)monomials(j, 0));
+ if (monomials.size2() > 1) p[j] *= pow_int(v, (int)monomials(j, 1));
+ if (monomials.size2() > 2) p[j] *= pow_int(w, (int)monomials(j, 2));
+ }
}
-}
-
-
+};
#endif
diff --git a/contrib/voro++/examples/degenerate/degenerate2.cc b/contrib/voro++/examples/degenerate/degenerate2.cc
index 6165d605d1..7ffe17a875 100644
--- a/contrib/voro++/examples/degenerate/degenerate2.cc
+++ b/contrib/voro++/examples/degenerate/degenerate2.cc
@@ -10,7 +10,7 @@ using namespace voro;
const double pi=3.1415926535897932384626433832795;
// The total number of points to create as degenerate vertices
-const int points=100;
+const int voro_points=100;
// The number of planes that will be cut around each point
const int n=64;
@@ -32,7 +32,7 @@ int main() {
v.init(-1,1,-1,1,-1,1);
// Plane cutting
- while(n<points) {
+ while(n<voro_points) {
// Choose a random point
x=2*rnd()-1;
--
GitLab