diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 65841628a6c3eac9b63465ac8acdbc81f2b0a0cd..5fcec7c8cb3b1460797887f9d72508f31ac0c2bb 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1209,10 +1209,9 @@ void GFace::moveToValidRange(SPoint2 &pt) const
}
}
-
-void GFace::addLayersOfQuads(int nLayers, GVertex *gv, double hmin, double ratio){
-
- SVector3 ez (0,0,1);
+void GFace::addLayersOfQuads(int nLayers, GVertex *gv, double hmin, double ratio)
+{
+ SVector3 ez (0, 0, 1);
std::list<GEdgeLoop>::iterator it = edgeLoops.begin();
FILE *f = fopen ("coucou.pos","w");
fprintf(f,"View \"\"{\n");
@@ -1231,21 +1230,23 @@ void GFace::addLayersOfQuads(int nLayers, GVertex *gv, double hmin, double ratio
GEdge *ge = it2->ge;
SPoint2 p[2];
if (it2->_sign == 1){
- for (int i=0;i<ge->lines.size();i++){
- reparamMeshEdgeOnFace ( ge->lines[i]->getVertex(0), ge->lines[i]->getVertex(1),this,p[0],p[1]);
- contour.push_back(std::make_pair(ge->lines[i]->getVertex(0),p[0]));
+ for (unsigned int i = 0; i < ge->lines.size(); i++){
+ reparamMeshEdgeOnFace(ge->lines[i]->getVertex(0), ge->lines[i]->getVertex(1),
+ this, p[0], p[1]);
+ contour.push_back(std::make_pair(ge->lines[i]->getVertex(0), p[0]));
}
}
else {
- for (int i=ge->lines.size()-1;i>=0;i--){
- reparamMeshEdgeOnFace ( ge->lines[i]->getVertex(0), ge->lines[i]->getVertex(1),this,p[0],p[1]);
+ for(int i = ge->lines.size() - 1; i >= 0; i--){
+ reparamMeshEdgeOnFace(ge->lines[i]->getVertex(0), ge->lines[i]->getVertex(1),
+ this,p[0],p[1]);
contour.push_back(std::make_pair(ge->lines[i]->getVertex(1),p[1]));
}
}
}
double hlayer = hmin;
- for (int j=0;j<nLayers;j++){
- for (int i=0;i<contour.size();i++){
+ for (int j = 0; j < nLayers; j++){
+ for (unsigned int i = 0; i < contour.size(); i++){
SPoint2 p0 = contour[(i+0) % contour.size()].second;
SPoint2 p1 = contour[(i+1) % contour.size()].second;
SPoint2 p2 = contour[(i+2) % contour.size()].second;
diff --git a/Geo/GeomMeshMatcher.cpp b/Geo/GeomMeshMatcher.cpp
index 1e8e56857482b05c6a454876937b2babea25683a..04b501d11d1cc4aaab391cf167eeb9521701644c 100644
--- a/Geo/GeomMeshMatcher.cpp
+++ b/Geo/GeomMeshMatcher.cpp
@@ -135,7 +135,7 @@ GeomMeshMatcher::matchVertices(GModel* m1, GModel *m2, bool& ok)
((GVertex*) best_candidate_ge)->setTag(((GVertex*) *entity1)->tag());
//m2->remove((GVertex*) best_candidate_ge);
//vertices.push_back((GVertex*) best_candidate_ge);
- for (int v = 0; v < ((GVertex*) best_candidate_ge)->getNumMeshVertices(); v++) {
+ for (unsigned int v = 0; v < ((GVertex*) best_candidate_ge)->getNumMeshVertices(); v++) {
((GVertex*) best_candidate_ge)->getMeshVertex(v)->setEntity((GVertex*) *entity1);
}
num_matched_vertices++;
@@ -225,7 +225,7 @@ GeomMeshMatcher::matchEdges(GModel* m1, GModel* m2,
choice->reverse();
}
- for (int v = 0; v < ((GEdge*) choice)->getNumMeshVertices(); v++) {
+ for (unsigned int v = 0; v < ((GEdge*) choice)->getNumMeshVertices(); v++) {
if (((GEdge*) choice)->getMeshVertex(v)->onWhat()->dim() > 0)
((GEdge*) choice)->getMeshVertex(v)->setEntity((GEdge*) *entity1);
}
@@ -295,7 +295,7 @@ GeomMeshMatcher:: matchFaces(GModel* m1, GModel* m2,
coresp_f->push_back(Pair<GFace*,GFace*>((GFace*) *entity1 ,
choice));
choice->setTag(((GFace*) *entity1)->tag());
- for (int v = 0; v < ((GFace*) choice)->getNumMeshVertices(); v++) {
+ for (unsigned int v = 0; v < ((GFace*) choice)->getNumMeshVertices(); v++) {
if(((GFace*) choice)->getMeshVertex(v)->onWhat()->dim() > 1)
((GFace*) choice)->getMeshVertex(v)->setEntity((GFace*) *entity1);
}
@@ -394,7 +394,7 @@ GeomMeshMatcher::matchRegions(GModel* m1, GModel* m2,
choice));
choice->setTag(((GRegion*) *entity1)->tag());
- for (int v = 0; v < ((GRegion*) choice)->getNumMeshVertices(); v++) {
+ for (unsigned int v = 0; v < ((GRegion*) choice)->getNumMeshVertices(); v++) {
if ( ((GRegion*) choice)->getMeshVertex(v)->onWhat()->dim() > 2)
((GRegion*) choice)->getMeshVertex(v)->setEntity((GRegion*) *entity1);
}
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 7bd56a56ad14997f4b05f21b99654ba06c16fb89..7b46f6ebdce65d3e989eff98956844cc57a02635 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -475,7 +475,8 @@ double MElement::integrate(double val[], int pOrder, int stride, int order)
return sum;
}
-static int getTriangleType (int order) {
+static int getTriangleType (int order)
+{
switch(order) {
case 0 : return MSH_TRI_1;
case 1 : return MSH_TRI_3;
@@ -488,10 +489,12 @@ static int getTriangleType (int order) {
case 8 : return MSH_TRI_45;
case 9 : return MSH_TRI_55;
case 10 : return MSH_TRI_66;
- default : Msg::Error("triangle order %i unknown", order);
+ default : Msg::Error("triangle order %i unknown", order); return 0;
}
}
-static int getQuadType (int order) {
+
+static int getQuadType (int order)
+{
switch(order) {
case 0 : return MSH_QUA_1;
case 1 : return MSH_QUA_4;
@@ -504,10 +507,12 @@ static int getQuadType (int order) {
case 8 : return MSH_QUA_81;
case 9 : return MSH_QUA_100;
case 10 : return MSH_QUA_121;
- default : Msg::Error("quad order %i unknown", order);
+ default : Msg::Error("quad order %i unknown", order); return 0;
}
}
-static int getLineType (int order) {
+
+static int getLineType (int order)
+{
switch(order) {
case 0 : return MSH_LIN_1;
case 1 : return MSH_LIN_2;
@@ -520,7 +525,7 @@ static int getLineType (int order) {
case 8 : return MSH_LIN_9;
case 9 : return MSH_LIN_10;
case 10 : return MSH_LIN_11;
- default : Msg::Error("line order %i unknown", order);
+ default : Msg::Error("line order %i unknown", order); return 0;
}
}
diff --git a/Graphics/Camera.h b/Graphics/Camera.h
index 090377f0f706670d6cede2ea456cf0df2d266f10..2c95f14d6b891d0b3a7d1f1d08be2f71b887a7f3 100644
--- a/Graphics/Camera.h
+++ b/Graphics/Camera.h
@@ -59,7 +59,7 @@ class Camera {
bool stereoEnable;
double Lc, eye_sep_ratio, closeness, ndfl, glFnear, glFfar, radians, wd2;
double glFleft,glFright,glFtop,glFbottom;
- Camera() : stereoEnable(false), on(false) {}
+ Camera() : on(false), stereoEnable(false) {}
~Camera(){}
void giveViewportDimension(const int& W,const int& H);
void lookAtCg();
diff --git a/Mesh/BoundaryLayers.cpp b/Mesh/BoundaryLayers.cpp
index 3c5e6dfdf9be3d3a62204a9328b4ebeddb4942c6..153b11ad0989af13f9032ad0c51b625153f62f6f 100644
--- a/Mesh/BoundaryLayers.cpp
+++ b/Mesh/BoundaryLayers.cpp
@@ -82,7 +82,7 @@ static void addExtrudeNormals(std::set<T*> &entities,
OctreePost *octree = 0;
#if defined(HAVE_POST)
if(view != -1){
- if(view >= 0 && view < PView::list.size()){
+ if(view >= 0 && view < (int)PView::list.size()){
octree = new OctreePost(PView::list[view]);
if(PView::list[view]->getData()->getNumVectors())
gouraud = false;
@@ -109,7 +109,7 @@ static void addExtrudeNormals(std::set<T*> &entities,
}
// enforce coherent normals at some points if necessary
- for(int i = 0; i < ExtrudeParams::normalsCoherence.size(); i++){
+ for(unsigned int i = 0; i < ExtrudeParams::normalsCoherence.size(); i++){
SPoint3 &p(ExtrudeParams::normalsCoherence[i]);
double n0[3], n1[3];
ExtrudeParams::normals[0]->get(p.x(), p.y(), p.z(), 3, n0);
diff --git a/Numeric/DivideAndConquer.cpp b/Numeric/DivideAndConquer.cpp
index b1de08dcf06404c7b754baad3bef56d0b64a5265..e732bcd9b0aafeb26d10c0127883328592ef33a6 100644
--- a/Numeric/DivideAndConquer.cpp
+++ b/Numeric/DivideAndConquer.cpp
@@ -866,8 +866,8 @@ void DocRecord::setPoints(fullMatrix<double> *p)
void DocRecord::initialize()
{
int i;
- for(i=0;i<numPoints;i++){
- points[i].flag = 0;
+ for(i = 0; i < numPoints; i++){
+ points[i].flag = 0;
}
}
@@ -888,21 +888,21 @@ void DocRecord::remove_all()
PointRecord* points2;
numPoints2 = 0;
for(i=0;i<numPoints;i++){
- if(points[i].flag==0){
+ if(points[i].flag == 0){
numPoints2++;
}
}
points2 = new PointRecord[numPoints2];
index = 0;
- for(i=0;i<numPoints;i++){
- if(points[i].flag==0){
- points2[index].where.h = points[i].where.h;
- points2[index].where.v = points[i].where.v;
- points2[index].data = points[i].data;
- points2[index].flag = points[i].flag;
- points2[index].identificator = points[i].identificator;
- index++;
- }
+ for(i = 0; i < numPoints; i++){
+ if(points[i].flag==0){
+ points2[index].where.h = points[i].where.h;
+ points2[index].where.v = points[i].where.v;
+ points2[index].data = points[i].data;
+ points2[index].flag = points[i].flag;
+ points2[index].identificator = points[i].identificator;
+ index++;
+ }
}
delete [] points;
points = points2;
diff --git a/Numeric/DivideAndConquer.h b/Numeric/DivideAndConquer.h
index b1ffec78c7965cc994fcdecbf5679a9b9bf1c21c..98c6b246221f77fe5631977743053c9d8a077d58 100644
--- a/Numeric/DivideAndConquer.h
+++ b/Numeric/DivideAndConquer.h
@@ -30,7 +30,7 @@ struct PointRecord {
void *data;
int flag; //0:to be kept, 1:to be removed
int identificator;
- PointRecord() : adjacent(0), data (0) {}
+ PointRecord() : adjacent(0), data (0), flag(0) {}
};
struct _CDLIST{
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index 48c164af035a93f0f7bee070ba6b3843b7f7eb6e..92412391cfb6c11ef38506df9176a98404f0e578 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -1,1391 +1,1380 @@
-// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to <gmsh@geuz.org>.
-
-#include "GmshDefines.h"
-#include "GmshMessage.h"
-#include "JacobianBasis.h"
-#include <vector>
-#include "polynomialBasis.h"
-
-// Bezier Exposants
-static fullMatrix<double> generate1DExposants(int order)
-{
- fullMatrix<double> exposants(order + 1, 1);
- exposants(0, 0) = 0;
- if (order > 0) {
- exposants(1, 0) = order;
- for (int i = 2; i < order + 1; i++) exposants(i, 0) = i - 1;
- }
-
- return exposants;
-}
-
-static fullMatrix<double> generateExposantsTriangle(int order)
-{
- int nbPoints = (order + 1) * (order + 2) / 2;
- fullMatrix<double> exposants(nbPoints, 2);
-
- exposants(0, 0) = 0;
- exposants(0, 1) = 0;
-
- if (order > 0) {
- exposants(1, 0) = order;
- exposants(1, 1) = 0;
- exposants(2, 0) = 0;
- exposants(2, 1) = order;
-
- if (order > 1) {
- int index = 3, ksi = 0, eta = 0;
-
- for (int i = 0; i < order - 1; i++, index++) {
- ksi++;
- exposants(index, 0) = ksi;
- exposants(index, 1) = eta;
- }
-
- ksi = order;
-
- for (int i = 0; i < order - 1; i++, index++) {
- ksi--;
- eta++;
- exposants(index, 0) = ksi;
- exposants(index, 1) = eta;
- }
-
- eta = order;
- ksi = 0;
-
- for (int i = 0; i < order - 1; i++, index++) {
- eta--;
- exposants(index, 0) = ksi;
- exposants(index, 1) = eta;
- }
-
- if (order > 2) {
- fullMatrix<double> inner = generateExposantsTriangle(order - 3);
- inner.add(1);
- exposants.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
- }
- }
- }
-
- return exposants;
-}
-static fullMatrix<double> generateExposantsQuad(int order)
-{
- int nbPoints = (order+1)*(order+1);
- fullMatrix<double> exposants(nbPoints, 2);
-
- exposants(0, 0) = 0;
- exposants(0, 1) = 0;
-
- if (order > 0) {
- exposants(1, 0) = order;
- exposants(1, 1) = 0;
- exposants(2, 0) = order;
- exposants(2, 1) = order;
- exposants(3, 0) = 0;
- exposants(3, 1) = order;
-
- if (order > 1) {
- int index = 4;
- const static int edges[4][2]={{0,1},{1,2},{2,3},{3,0}};
- for (int iedge=0; iedge<4; iedge++) {
- int p0 = edges[iedge][0];
- int p1 = edges[iedge][1];
- for (int i = 1; i < order; i++, index++) {
- exposants(index, 0) = exposants(p0, 0) + i*(exposants(p1,0)-exposants(p0,0))/order;
- exposants(index, 1) = exposants(p0, 1) + i*(exposants(p1,1)-exposants(p0,1))/order;
- }
- }
- if (order > 2) {
- fullMatrix<double> inner = generateExposantsQuad(order - 2);
- inner.add(1);
- exposants.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
- }
- }
- }
-
- return exposants;
-}
-static int nbdoftriangle(int order) { return (order + 1) * (order + 2) / 2; }
-
-static void nodepositionface0(int order, double *u, double *v, double *w)
-{ // uv surface - orientation v0-v2-v1
- int ndofT = nbdoftriangle(order);
- if (order == 0) { u[0] = 0.; v[0] = 0.; w[0] = 0.; return; }
-
- u[0]= 0.; v[0]= 0.; w[0] = 0.;
- u[1]= 0.; v[1]= order; w[1] = 0.;
- u[2]= order; v[2]= 0.; w[2] = 0.;
-
- // 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] = k + 1;
- v[3 + order - 1 + k] = order - 1 - k ;
- w[3 + order - 1 + k] = 0.;
-
- u[3 + 2 * (order - 1) + k] = order - 1 - k;
- v[3 + 2 * (order - 1) + k] = 0.;
- 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);
- }
- }
- for (int k = 0; k < ndofT; k++){
- u[k] = u[k] / order;
- v[k] = v[k] / order;
- w[k] = w[k] / order;
- }
-}
-
-static void nodepositionface1(int order, double *u, double *v, double *w)
-{ // uw surface - orientation v0-v1-v3
- int ndofT = nbdoftriangle(order);
- if (order == 0) { u[0] = 0.; v[0] = 0.; w[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;
- }
- 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;
- }
-}
-
-static void nodepositionface2(int order, double *u, double *v, double *w)
-{ // vw surface - orientation v0-v3-v2
- int ndofT = nbdoftriangle(order);
- if (order == 0) { u[0] = 0.; v[0] = 0.; return; }
-
- u[0]= 0.; v[0]= 0.; w[0] = 0.;
- u[1]= 0.; v[1]= 0.; w[1] = order;
- u[2]= 0.; v[2]= order; w[2] = 0.;
- // edges
- for (int k = 0; k < (order - 1); k++){
-
- u[3 + k] = 0.;
- v[3 + k] = 0.;
- w[3 + k] = k + 1;
-
- u[3 + order - 1 + k] = 0.;
- v[3 + order - 1 + k] = k + 1;
- w[3 + order - 1 + k] = order - 1 - k;
-
- 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);
- 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);
- 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;
- }
-}
-
-static void nodepositionface3(int order, double *u, double *v, double *w)
-{ // uvw surface - orientation v3-v1-v2
- int ndofT = nbdoftriangle(order);
- if (order == 0) { u[0] = 0.; v[0] = 0.; w[0] = 0.; return; }
-
- u[0]= 0.; v[0]= 0.; w[0] = order;
- 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] = order - 1 - k;
-
- 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] = k + 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) + 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;
- }
-}
-
-static fullMatrix<double> generateExposantsTetrahedron(int order)
-{
- int nbPoints = (order + 1) * (order + 2) * (order + 3) / 6;
-
- fullMatrix<double> exposants(nbPoints, 3);
-
- exposants(0, 0) = 0;
- exposants(0, 1) = 0;
- exposants(0, 2) = 0;
-
- if (order > 0) {
- exposants(1, 0) = order;
- exposants(1, 1) = 0;
- exposants(1, 2) = 0;
-
- exposants(2, 0) = 0;
- exposants(2, 1) = order;
- exposants(2, 2) = 0;
-
- exposants(3, 0) = 0;
- exposants(3, 1) = 0;
- exposants(3, 2) = order;
-
- // edges e5 and e6 switched in original version, opposite direction
- // the template has been defined in table edges_tetra and faces_tetra (MElement.h)
-
- if (order > 1) {
- for (int k = 0; k < (order - 1); k++) {
- exposants(4 + k, 0) = k + 1;
- exposants(4 + order - 1 + k, 0) = order - 1 - k;
- exposants(4 + 2 * (order - 1) + k, 0) = 0;
- exposants(4 + 3 * (order - 1) + k, 0) = 0;
- // exposants(4 + 4 * (order - 1) + k, 0) = order - 1 - k;
- // exposants(4 + 5 * (order - 1) + k, 0) = 0.;
- exposants(4 + 4 * (order - 1) + k, 0) = 0;
- exposants(4 + 5 * (order - 1) + k, 0) = k+1;
-
- exposants(4 + k, 1) = 0;
- exposants(4 + order - 1 + k, 1) = k + 1;
- exposants(4 + 2 * (order - 1) + k, 1) = order - 1 - k;
- exposants(4 + 3 * (order - 1) + k, 1) = 0;
- // exposants(4 + 4 * (order - 1) + k, 1) = 0.;
- // exposants(4 + 5 * (order - 1) + k, 1) = order - 1 - k;
- exposants(4 + 4 * (order - 1) + k, 1) = k+1;
- exposants(4 + 5 * (order - 1) + k, 1) = 0;
-
- exposants(4 + k, 2) = 0;
- exposants(4 + order - 1 + k, 2) = 0;
- exposants(4 + 2 * (order - 1) + k, 2) = 0;
- exposants(4 + 3 * (order - 1) + k, 2) = order - 1 - k;
- exposants(4 + 4 * (order - 1) + k, 2) = order - 1 - k;
- exposants(4 + 5 * (order - 1) + k, 2) = order - 1 - k;
- }
-
- if (order > 2) {
- 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);
-
- // u-v plane
-
- for (int i = 0; i < nbdofface; i++){
- exposants(ns + i, 0) = u[i] * (order - 3) + 1;
- exposants(ns + i, 1) = v[i] * (order - 3) + 1;
- exposants(ns + i, 2) = w[i] * (order - 3);
- }
-
- ns = ns + nbdofface;
-
- // u-w plane
-
- nodepositionface1(order - 3, u, v, w);
-
- for (int i=0; i < nbdofface; i++){
- exposants(ns + i, 0) = u[i] * (order - 3) + 1;
- exposants(ns + i, 1) = v[i] * (order - 3) ;
- exposants(ns + i, 2) = w[i] * (order - 3) + 1;
- }
-
- // v-w plane
-
- ns = ns + nbdofface;
-
- nodepositionface2(order - 3, u, v, w);
-
- for (int i = 0; i < nbdofface; i++){
- exposants(ns + i, 0) = u[i] * (order - 3);
- exposants(ns + i, 1) = v[i] * (order - 3) + 1;
- exposants(ns + i, 2) = w[i] * (order - 3) + 1;
- }
-
- // u-v-w plane
-
- ns = ns + nbdofface;
-
- nodepositionface3(order - 3, u, v, w);
-
- for (int i = 0; i < nbdofface; i++){
- exposants(ns + i, 0) = u[i] * (order - 3) + 1;
- exposants(ns + i, 1) = v[i] * (order - 3) + 1;
- exposants(ns + i, 2) = w[i] * (order - 3) + 1;
- }
-
- ns = ns + nbdofface;
-
- delete [] u;
- delete [] v;
- delete [] w;
-
- if (order > 3) {
-
- fullMatrix<double> interior = generateExposantsTetrahedron(order - 4);
- for (int k = 0; k < interior.size1(); k++) {
- exposants(ns + k, 0) = 1 + interior(k, 0);
- exposants(ns + k, 1) = 1 + interior(k, 1);
- exposants(ns + k, 2) = 1 + interior(k, 2);
- }
- }
- }
- }
- }
-
- return exposants;
-}
-
-static fullMatrix<double> generateExposantsPrism(int order)
-{
-// int nbMonomials = (order + 1) * (order + 1) * (order + 2) / 2;
-//
-// fullMatrix<double> monomials(nbMonomials, 3);
-// int index = 0;
-// fullMatrix<double> lineMonoms = generate1DMonomials(order);
-// fullMatrix<double> triMonoms = generateExposantsTriangle(order);
-// // store monomials in right order
-// for (int currentOrder = 0; currentOrder <= order; currentOrder++) {
-// int orderT = currentOrder, orderL = currentOrder;
-// for (orderL = 0; orderL < currentOrder; orderL++) {
-// // do all permutations of monoms for orderL, orderT
-// int iL = orderL;
-// for (int iT = (orderT)*(orderT+1)/2; iT < (orderT+1)*(orderT+2)/2 ;iT++) {
-// monomials(index,0) = triMonoms(iT,0);
-// monomials(index,1) = triMonoms(iT,1);
-// monomials(index,2) = lineMonoms(iL,0);
-// index ++;
-// }
-// }
-// orderL = currentOrder;
-// for (orderT = 0; orderT <= currentOrder; orderT++) {
-// int iL = orderL;
-// for (int iT = (orderT)*(orderT+1)/2; iT < (orderT+1)*(orderT+2)/2 ;iT++) {
-// monomials(index,0) = triMonoms(iT,0);
-// monomials(index,1) = triMonoms(iT,1);
-// monomials(index,2) = lineMonoms(iL,0);
-// index ++;
-// }
-// }
-// }
-//// monomials.print("Pri monoms");
-// return monomials;
-
- //const double prism18Pts[18][3] = {
- // {0, 0, -1}, // 0
- // {1, 0, -1}, // 1
- // {0, 1, -1}, // 2
- // {0, 0, 1}, // 3
- // {1, 0, 1}, // 4
- // {0, 1, 1}, // 5
- // {0.5, 0, -1}, // 6
- // {0, 0.5, -1}, // 7
- // {0, 0, 0}, // 8
- // {0.5, 0.5, -1}, // 9
- // {1, 0, 0}, // 10
- // {0, 1, 0}, // 11
- // {0.5, 0, 1}, // 12
- // {0, 0.5, 1}, // 13
- // {0.5, 0.5, 1}, // 14
- // {0.5, 0, 0}, // 15
- // {0, 0.5, 0}, // 16
- // {0.5, 0.5, 0}, // 17
- // };
-
- int nbPoints = (order + 1)*(order + 1)*(order + 2)/2;
- fullMatrix<double> exposants(nbPoints, 3);
-
- int index = 0;
- fullMatrix<double> triExp = generateExposantsTriangle(order);
- fullMatrix<double> lineExp = generate1DExposants(order);
-
- /* if (order == 2)
- for (int i =0; i<18; i++)
- for (int j=0; j<3;j++)
- exposants(i,j) = prism18Pts[i][j];
- else*/
- {
- for (int i = 0; i < 3; i++) {
- exposants(index,0) = triExp(i,0);
- exposants(index,1) = triExp(i,1);
- exposants(index,2) = lineExp(0,0);
- index ++;
- exposants(index,0) = triExp(i,0);
- exposants(index,1) = triExp(i,1);
- exposants(index,2) = lineExp(1,0);
- index ++;
- }
- for (int i = 3; i < triExp.size1(); i++) {
- exposants(index,0) = triExp(i,0);
- exposants(index,1) = triExp(i,1);
- exposants(index,2) = lineExp(0,0);
- index ++;
- exposants(index,0) = triExp(i,0);
- exposants(index,1) = triExp(i,1);
- exposants(index,2) = lineExp(1,0);
- index ++;
- }
- for (int j = 2; j <lineExp.size1() ; j++) {
- for (int i = 0; i < triExp.size1(); i++) {
- exposants(index,0) = triExp(i,0);
- exposants(index,1) = triExp(i,1);
- exposants(index,2) = lineExp(j,0);
- index ++;
- }
- }
- }
-
- return exposants;
-}
-
-// Sampling Points
-static fullMatrix<double> generate1DPoints(int order)
-{
- fullMatrix<double> line(order + 1, 1);
- line(0,0) = 0;
- if (order > 0) {
- line(0, 0) = 0.;
- line(1, 0) = 1.;
- double dd = 1. / order;
- for (int i = 2; i < order + 1; i++) line(i, 0) = dd * (i - 1);
- }
-
- return line;
-}
-
-static fullMatrix<double> 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);
-
- fullMatrix<double> point(nbPoints, 3);
-
- double overOrder = (order == 0 ? 1. : 1. / order);
-
- 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(2, 0) = 0.;
- point(2, 1) = order;
- point(2, 2) = 0.;
-
- point(3, 0) = 0.;
- point(3, 1) = 0.;
- point(3, 2) = order;
-
- // edges e5 and e6 switched in original version, opposite direction
- // the template has been defined in table edges_tetra and faces_tetra (MElement.h)
-
- 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.;
- point(4 + 4 * (order - 1) + k, 0) = 0.;
- point(4 + 5 * (order - 1) + k, 0) = k+1;
-
- 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 + 4 * (order - 1) + k, 1) = k+1;
- point(4 + 5 * (order - 1) + k, 1) = 0.;
-
- 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) = order - 1 - k;
- point(4 + 4 * (order - 1) + k, 2) = order - 1 - k;
- point(4 + 5 * (order - 1) + k, 2) = order - 1 - k;
- }
-
- if (order > 2) {
- 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);
-
- // u-v plane
-
- 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;
-
- // u-w plane
-
- 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.;
- }
-
- // v-w plane
-
- ns = ns + nbdofface;
-
- nodepositionface2(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.;
- }
-
- // u-v-w plane
-
- ns = ns + nbdofface;
-
- nodepositionface3(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.;
- }
-
- ns = ns + nbdofface;
-
- delete [] u;
- delete [] v;
- delete [] w;
-
- if (!serendip && order > 3) {
-
- fullMatrix<double> interior = gmshGeneratePointsTetrahedron(order - 4, false);
- for (int 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);
- }
- }
- }
- }
- }
-
- point.scale(overOrder);
- return point;
-}
-
-static fullMatrix<double> gmshGeneratePointsTriangle(int order, bool serendip)
-{
- int nbPoints = serendip ? 3 * order : (order + 1) * (order + 2) / 2;
- fullMatrix<double> point(nbPoints, 2);
-
- point(0, 0) = 0;
- point(0, 1) = 0;
-
- if (order > 0) {
- double dd = 1. / order;
-
- point(1, 0) = 1;
- point(1, 1) = 0;
- point(2, 0) = 0;
- point(2, 1) = 1;
-
- int index = 3;
-
- if (order > 1) {
-
- double ksi = 0;
- double eta = 0;
-
- for (int i = 0; i < order - 1; i++, index++) {
- ksi += dd;
- point(index, 0) = ksi;
- point(index, 1) = eta;
- }
-
- ksi = 1.;
-
- for (int i = 0; i < order - 1; i++, index++) {
- ksi -= dd;
- eta += dd;
- point(index, 0) = ksi;
- point(index, 1) = eta;
- }
-
- eta = 1.;
- ksi = 0.;
-
- for (int i = 0; i < order - 1; i++, index++) {
- eta -= dd;
- point(index, 0) = ksi;
- point(index, 1) = eta;
- }
-
- if (order > 2 && !serendip) {
- fullMatrix<double> inner = gmshGeneratePointsTriangle(order - 3, serendip);
- inner.scale(1. - 3. * dd);
- inner.add(dd);
- point.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
- }
- }
- }
- return point;
-}
-
-static fullMatrix<double> generatePointsPrism(int order, bool serendip)
-{
- const double prism18Pts[18][3] = {
- {0, 0, 0}, // 0
- {1, 0, 0}, // 1
- {0, 1, 0}, // 2
- {0, 0, 1}, // 3
- {1, 0, 1}, // 4
- {0, 1, 1}, // 5
- {.5, 0, 0}, // 6
- {0, .5, 0}, // 7
- {0, 0, .5}, // 8
- {.5, .5, 0}, // 9
- {1, 0, .5}, // 10
- {0, 1, .5}, // 11
- {.5, 0, 1}, // 12
- {0, .5, 1}, // 13
- {.5, .5, 1}, // 14
- {.5, 0, .5}, // 15
- {0, .5, .5}, // 16
- {.5, .5, .5}, // 17
- };
-
- int nbPoints = (order + 1)*(order + 1)*(order + 2)/2;
- fullMatrix<double> point(nbPoints, 3);
-
- int index = 0;
-
- if (order == 2)
- for (int i =0; i<18; i++)
- for (int j=0; j<3;j++)
- point(i,j) = prism18Pts[i][j];
- else {
- fullMatrix<double> triPoints = gmshGeneratePointsTriangle(order,false);
- fullMatrix<double> linePoints = generate1DPoints(order);
- for (int j = 0; j <linePoints.size1() ; j++) {
- for (int i = 0; i < triPoints.size1(); i++) {
- point(index,0) = triPoints(i,0);
- point(index,1) = triPoints(i,1);
- point(index,2) = linePoints(j,0);
- index ++;
- }
- }
- }
- return point;
-}
-
-static fullMatrix<double> generatePointsQuad(int order, bool serendip)
-{
- int nbPoints = serendip ? order*4 : (order+1)*(order+1);
- fullMatrix<double> point(nbPoints, 2);
-
- if (order > 0) {
- point(0, 0) = 0;
- point(0, 1) = 0;
- point(1, 0) = 1;
- point(1, 1) = 0;
- point(2, 0) = 1;
- point(2, 1) = 1;
- point(3, 0) = 0;
- point(3, 1) = 1;
-
- if (order > 1) {
- int index = 4;
- const static int edges[4][2]={{0,1},{1,2},{2,3},{3,0}};
- for (int iedge=0; iedge<4; iedge++) {
- int p0 = edges[iedge][0];
- int p1 = edges[iedge][1];
- for (int i = 1; i < order; i++, index++) {
- point(index, 0) = point(p0, 0) + i*(point(p1,0)-point(p0,0))/order;
- point(index, 1) = point(p0, 1) + i*(point(p1,1)-point(p0,1))/order;
- }
- }
- if (order > 2 && !serendip) {
- fullMatrix<double> inner = generatePointsQuad(order - 2, false);
- inner.scale(1. - 2./order);
- point.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
- }
- }
- }
- else {
- point(0, 0) = 0;
- point(0, 1) = 0;
- }
- return point;
-}
-
-// Sub Control Points
-static std::vector< fullMatrix<double> > generateSubPointsLine(int order)
-{
- int nbPts = (order + 1);
-
- std::vector< fullMatrix<double> > subPoints(2);
- fullMatrix<double> prox;
-
- subPoints[0] = generate1DPoints(order);
- subPoints[0].scale(.5);
-
- subPoints[1].copy(subPoints[0]);
- subPoints[1].add(.5);
-
- return subPoints;
-}
-
-static std::vector< fullMatrix<double> > generateSubPointsTriangle(int order)
-{
- int nbPts = (order + 1) * (order + 2) / 2;
-
- std::vector< fullMatrix<double> > subPoints(4);
- fullMatrix<double> prox;
-
- subPoints[0] = gmshGeneratePointsTriangle(order, false);
- subPoints[0].scale(.5);
-
- subPoints[1].copy(subPoints[0]);
- prox.setAsProxy(subPoints[1], 0, 1);
- prox.add(.5);
-
- subPoints[2].copy(subPoints[0]);
- prox.setAsProxy(subPoints[2], 1, 1);
- prox.add(.5);
-
- subPoints[3].copy(subPoints[0]);
- subPoints[3].scale(-1.);
- subPoints[3].add(.5);
-
- return subPoints;
-}
-
-static std::vector< fullMatrix<double> > generateSubPointsTetrahedron(int order)
-{
- int nbPts = (order + 1) * (order + 2) * (order + 3) / 6;
-
- std::vector< fullMatrix<double> > subPoints(8);
- fullMatrix<double> prox1;
- fullMatrix<double> prox2;
-
- subPoints[0] = gmshGeneratePointsTetrahedron(order, false);
- subPoints[0].scale(.5);
-
- subPoints[1].copy(subPoints[0]);
- prox1.setAsProxy(subPoints[1], 0, 1);
- prox1.add(.5);
-
- subPoints[2].copy(subPoints[0]);
- prox1.setAsProxy(subPoints[2], 1, 1);
- prox1.add(.5);
-
- subPoints[3].copy(subPoints[0]);
- prox1.setAsProxy(subPoints[3], 2, 1);
- prox1.add(.5);
-
- // u := .5-u-w
- // v := .5-v-w
- // w := w
- subPoints[4].copy(subPoints[0]);
- prox1.setAsProxy(subPoints[4], 0, 2);
- prox1.scale(-1.);
- prox1.add(.5);
- prox1.setAsProxy(subPoints[4], 0, 1);
- prox2.setAsProxy(subPoints[4], 2, 1);
- prox1.add(prox2, -1.);
- prox1.setAsProxy(subPoints[4], 1, 1);
- prox1.add(prox2, -1.);
-
- // u := u
- // v := .5-v
- // w := w+v
- subPoints[5].copy(subPoints[0]);
- prox1.setAsProxy(subPoints[5], 2, 1);
- prox2.setAsProxy(subPoints[5], 1, 1);
- prox1.add(prox2);
- prox2.scale(-1.);
- prox2.add(.5);
-
- // u := .5-u
- // v := v
- // w := w+u
- subPoints[6].copy(subPoints[0]);
- prox1.setAsProxy(subPoints[6], 2, 1);
- prox2.setAsProxy(subPoints[6], 0, 1);
- prox1.add(prox2);
- prox2.scale(-1.);
- prox2.add(.5);
-
- // u := u+w
- // v := v+w
- // w := .5-w
- subPoints[7].copy(subPoints[0]);
- prox1.setAsProxy(subPoints[7], 0, 1);
- prox2.setAsProxy(subPoints[7], 2, 1);
- prox1.add(prox2);
- prox1.setAsProxy(subPoints[7], 1, 1);
- prox1.add(prox2);
- prox2.scale(-1.);
- prox2.add(.5);
-
-
- return subPoints;
-}
-
-static std::vector< fullMatrix<double> > generateSubPointsQuad(int order)
-{
- int nbPts = (order + 1) * (order + 1);
-
- std::vector< fullMatrix<double> > subPoints(4);
- fullMatrix<double> prox;
-
- subPoints[0] = generatePointsQuad(order, false);
- subPoints[0].scale(.5);
-
- subPoints[1].copy(subPoints[0]);
- prox.setAsProxy(subPoints[1], 0, 1);
- prox.add(.5);
-
- subPoints[2].copy(subPoints[0]);
- prox.setAsProxy(subPoints[2], 1, 1);
- prox.add(.5);
-
- subPoints[3].copy(subPoints[1]);
- prox.setAsProxy(subPoints[3], 1, 1);
- prox.add(.5);
-
- return subPoints;
-}
-
-static std::vector< fullMatrix<double> > generateSubPointsPrism(int order)
-{
- int nbPts = (order + 1) * (order + 1) * (order + 2) / 2;
-
- std::vector< fullMatrix<double> > subPoints(8);
- fullMatrix<double> prox;
-
- subPoints[0] = generatePointsPrism(order, false);
- subPoints[0].scale(.5);
-
- subPoints[1].copy(subPoints[0]);
- prox.setAsProxy(subPoints[1], 0, 1);
- prox.add(.5);
-
- subPoints[2].copy(subPoints[0]);
- prox.setAsProxy(subPoints[2], 1, 1);
- prox.add(.5);
-
- subPoints[3].copy(subPoints[0]);
- prox.setAsProxy(subPoints[3], 0, 2);
- prox.scale(-1.);
- prox.add(.5);
-
- subPoints[4].copy(subPoints[0]);
- prox.setAsProxy(subPoints[4], 2, 1);
- prox.add(.5);
-
- subPoints[5].copy(subPoints[1]);
- prox.setAsProxy(subPoints[5], 2, 1);
- prox.add(.5);
-
- subPoints[6].copy(subPoints[2]);
- prox.setAsProxy(subPoints[6], 2, 1);
- prox.add(.5);
-
- subPoints[7].copy(subPoints[3]);
- prox.setAsProxy(subPoints[7], 2, 1);
- prox.add(.5);
-
- return subPoints;
-}
-
-// Matrices generation
-static int nChoosek(int n, int k)
-{
- if (n < k || k < 0) {
- Msg::Error("Wrong argument for combination.");
- return 1;
- }
-
- if (k > n/2) k = n-k;
- if (k == 1)
- return n;
- if (k == 0)
- return 1;
-
- int c = 1;
- for (int i = 1; i <= k; i++, n--) (c *= n) /= i;
- return c;
-}
-
-
-static fullMatrix<double> generateLag2BezMatrix
- (const fullMatrix<double> &exposant, const fullMatrix<double> &point,
- int order, int dimSimplex, bool invert = true)
-{
-
- if(exposant.size1() != point.size1() || exposant.size2() != point.size2()){
- Msg::Fatal("Wrong sizes for Bezier coefficients generation %d %d -- %d %d",
- exposant.size1(),point.size1(),
- exposant.size2(),point.size2());
- return fullMatrix<double>(1, 1);
- }
-
- int ndofs = exposant.size1();
- int dim = exposant.size2();
-
- fullMatrix<double> Vandermonde(ndofs, ndofs);
- for (int i = 0; i < ndofs; i++) {
- for (int j = 0; j < ndofs; j++) {
- double dd = 1.;
-
- double pointCompl = 1.;
- int exposantCompl = order;
- for (int k = 0; k < dimSimplex; k++) {
- dd *= nChoosek(exposantCompl, (int) exposant(i, k))
- * pow(point(j, k), exposant(i, k));
- pointCompl -= point(j, k);
- exposantCompl -= (int) exposant(i, k);
- }
- dd *= pow(pointCompl, exposantCompl);
-
- for (int k = dimSimplex; k < dim; k++)
- dd *= nChoosek(order, (int) exposant(i, k))
- * pow(point(j, k), exposant(i, k))
- * pow(1. - point(j, k), order - exposant(i, k));
-
- Vandermonde(j, i) = dd;
- }
- }
-
- if (!invert) return Vandermonde;
-
- fullMatrix<double> coefficient(ndofs, ndofs);
- Vandermonde.invert(coefficient);
- return coefficient;
-}
-
-
-
-static fullMatrix<double> generateDivisor
- (const fullMatrix<double> &exposants, const std::vector< fullMatrix<double> > &subPoints,
- const fullMatrix<double> &lag2Bez, int order, int dimSimplex)
-{
- if (exposants.size1() != lag2Bez.size1() || exposants.size1() != lag2Bez.size2()){
- Msg::Fatal("Wrong sizes for Bezier Divisor %d %d -- %d %d",
- exposants.size1(), lag2Bez.size1(),
- exposants.size1(), lag2Bez.size2());
- return fullMatrix<double>(1, 1);
- }
-
- int dim = exposants.size2();
- int nbPts = lag2Bez.size1();
- int nbSubPts = nbPts * subPoints.size();
-
- fullMatrix<double> intermediate2(nbPts, nbPts);
- fullMatrix<double> divisor(nbSubPts, nbPts);
-
- for (unsigned int i = 0; i < subPoints.size(); i++) {
- fullMatrix<double> intermediate1 =
- generateLag2BezMatrix(exposants, subPoints[i], order, dimSimplex, false);
- lag2Bez.mult(intermediate1, intermediate2);
- divisor.copy(intermediate2, 0, nbPts, 0, nbPts, i*nbPts, 0);
- }
- return divisor;
-}
-
-
-
-static void generateGradShapes(JacobianBasis &jfs, const fullMatrix<double> &points
- , const fullMatrix<double> &monomials, const fullMatrix<double> &coefficients)
-{
-
- /*{
- Msg::Info("Printing vector jacobian':");
- int ni = points.size1();
- int nj = points.size2();
- Msg::Info(" ");
- for(int I = 0; I < ni; I++){
- Msg::Info("%lf - %lf", points(I, 0), points(I, 1));
- }
- Msg::Info(" ");
- }
- {
- Msg::Info("Printing vector jacobian':");
- int ni = monomials.size1();
- int nj = monomials.size2();
- Msg::Info(" ");
- for(int I = 0; I < ni; I++){
- Msg::Info("%lf - %lf", monomials(I, 0), monomials(I, 1));
- }
- Msg::Info(" ");
- }
- {
- Msg::Info("Printing vector jacobian':");
- int ni = coefficients.size1();
- int nj = coefficients.size2();
- Msg::Info(" ");
- for(int I = 0; I < ni; I++){
- Msg::Info(" ");
- for(int J = 0; J < nj; J++){
- Msg::Info("%lf", coefficients(J, I));
- }
- }
- Msg::Info(" ");
- }*/
-
- int nbPts = points.size1();
- int nbDofs = monomials.size1();
- int dim = points.size2();
-
- switch (dim) {
- case 3 :
- jfs.gradShapeMatZ.resize(nbPts, nbDofs, true);
- case 2 :
- jfs.gradShapeMatY.resize(nbPts, nbDofs, true);
- case 1 :
- jfs.gradShapeMatX.resize(nbPts, nbDofs, true);
- break;
- default :
- return;
- }
-
- double dx, dy, dz;
-
- switch (dim) {
- case 3 :
- for (int i = 0; i < nbDofs; i++) {
- for (int k = 0; k < nbPts; k++) {
-
- if ((int) monomials(i, 0) > 0) {
- dx = pow( points(k, 0), monomials(i, 0)-1 ) * monomials(i, 0)
- * pow( points(k, 1), monomials(i, 1) )
- * pow( points(k, 2), monomials(i, 2) );
- for (int j = 0; j < nbDofs; j++)
- jfs.gradShapeMatX(k, j) += coefficients(j, i) * dx;
- }
- if ((int) monomials(i, 1) > 0.) {
- dy = pow( points(k, 0), monomials(i, 0) )
- * pow( points(k, 1), monomials(i, 1)-1 ) * monomials(i, 1)
- * pow( points(k, 2), monomials(i, 2) );
- for (int j = 0; j < nbDofs; j++)
- jfs.gradShapeMatY(k, j) += coefficients(j, i) * dy;
- }
- if ((int) monomials(i, 2) > 0.) {
- dz = pow( points(k, 0), monomials(i, 0) )
- * pow( points(k, 1), monomials(i, 1) )
- * pow( points(k, 2), monomials(i, 2)-1 ) * monomials(i, 2);
- for (int j = 0; j < nbDofs; j++)
- jfs.gradShapeMatZ(k, j) += coefficients(j, i) * dz;
- }
- }
- }
- return;
-
- case 2 :
- for (int i = 0; i < nbDofs; i++) {
- for (int k = 0; k < nbPts; k++) {
-
- if ((int) monomials(i, 0) > 0) {
- dx = pow( points(k, 0), (int) monomials(i, 0)-1 ) * monomials(i, 0)
- * pow( points(k, 1), (int) monomials(i, 1) );
- for (int j = 0; j < nbDofs; j++)
- jfs.gradShapeMatX(k, j) += coefficients(j, i) * dx;
- }
- if ((int) monomials(i, 1) > 0) {
- dy = pow( points(k, 0), (int) monomials(i, 0) )
- * pow( points(k, 1), (int) monomials(i, 1)-1 ) * monomials(i, 1);
- for (int j = 0; j < nbDofs; j++)
- jfs.gradShapeMatY(k, j) += coefficients(j, i) * dy;
- }
- }
- }
- return;
-
- case 1 :
- for (int i = 0; i < nbDofs; i++) {
- for (int k = 0; k < nbPts; k++) {
-
- if ((int) monomials(i, 0) > 0) {
- dx = pow( points(k, 0), (int) monomials(i, 0)-1 ) * monomials(i, 0);
- for (int j = 0; j < nbDofs; j++)
- jfs.gradShapeMatX(k, j) += coefficients(j, i) * dx;
- }
- }
- }
- }
- return;
-}
-
-std::map<int, JacobianBasis> JacobianBases::fs;
-
-const JacobianBasis *JacobianBases::find(int tag)
-{
- std::map<int, JacobianBasis>::const_iterator it = fs.find(tag);
- if (it != fs.end()) return &it->second;
-
- JacobianBasis B;
- B.numLagPts = -1;
-
- int order;
-
-
- if (tag == MSH_PNT) {
- B.numLagPts = 1;
- B.exposants = generate1DExposants(0);
- B.points = generate1DPoints(0);
- B.matrixLag2Bez = generateLag2BezMatrix(B.exposants, B.points, 0, 0);
- /*std::vector< fullMatrix<double> > subPoints = generateSubPointsLine(0);
- B.divisor = generateDivisor(B.exposants, subPoints, B.matrixLag2Bez, 0, 0);
- const polynomialBasis *F = polynomialBases::find(tag);
- generateGradShapes(B, B.points, F->monomials, F->coefficients);*/
- goto end;
- }
-
-
- switch (tag) {
- case MSH_LIN_2: order = 1; break;
- case MSH_LIN_3: order = 2; break;
- case MSH_LIN_4: order = 3; break;
- case MSH_LIN_5: order = 4; break;
- case MSH_LIN_6: order = 5; break;
- case MSH_LIN_7: order = 6; break;
- case MSH_LIN_8: order = 7; break;
- case MSH_LIN_9: order = 8; break;
- case MSH_LIN_10: order = 9; break;
- case MSH_LIN_11: order = 10; break;
- default: order = -1; break;
- }
- if (order >= 0) {
- int o = order - 1;
- B.numLagPts = 2;
- B.exposants = generate1DExposants(o);
- B.points = generate1DPoints(o);
- B.matrixLag2Bez = generateLag2BezMatrix(B.exposants, B.points, o, 0);
- std::vector< fullMatrix<double> > subPoints = generateSubPointsLine(0);
- B.divisor = generateDivisor(B.exposants, subPoints, B.matrixLag2Bez, o, 0);
- const polynomialBasis *F = polynomialBases::find(tag);
- generateGradShapes(B, B.points, F->monomials, F->coefficients);
- goto end;
- }
-
-
-
- switch (tag) {
- case MSH_TRI_3 : order = 1; break;
- case MSH_TRI_6 : order = 2; break;
- case MSH_TRI_9 :
- case MSH_TRI_10 : order = 3; break;
- case MSH_TRI_12 :
- case MSH_TRI_15 : order = 4; break;
- case MSH_TRI_15I :
- case MSH_TRI_21 : order = 5; break;
- case MSH_TRI_18 :
- case MSH_TRI_28 : order = 6; break;
- case MSH_TRI_21I :
- case MSH_TRI_36 : order = 7; break;
- case MSH_TRI_24 :
- case MSH_TRI_45 : order = 8; break;
- case MSH_TRI_27 :
- case MSH_TRI_55 : order = 9; break;
- case MSH_TRI_30 :
- case MSH_TRI_66 : order = 10; break;
- default: order = -1; break;
- }
- if (order >= 0) {
- int o = 2 * (order - 1);
- B.numLagPts = 3;
- B.exposants = generateExposantsTriangle(o);
- B.points = gmshGeneratePointsTriangle(o,false);
- B.matrixLag2Bez = generateLag2BezMatrix(B.exposants, B.points, o, 2);
- std::vector< fullMatrix<double> > subPoints = generateSubPointsTriangle(o);
- B.divisor = generateDivisor(B.exposants, subPoints, B.matrixLag2Bez, o, 2);
- const polynomialBasis *F = polynomialBases::find(tag);
- generateGradShapes(B, B.points, F->monomials, F->coefficients);
- goto end;
- }
-
-
- switch (tag) {
- case MSH_TET_4 : order = 1; break;
- case MSH_TET_10 : order = 2; break;
- case MSH_TET_20 : order = 3; break;
- case MSH_TET_34 :
- case MSH_TET_35 : order = 4; break;
- case MSH_TET_52 :
- case MSH_TET_56 : order = 5; break;
- case MSH_TET_84 : order = 6; break;
- case MSH_TET_120 : order = 7; break;
- case MSH_TET_165 : order = 8; break;
- case MSH_TET_220 : order = 9; break;
- case MSH_TET_286 : order = 10; break;
- default: order = -1; break;
- }
- if (order >= 0) {
- int o = 3 * (order - 1);
- B.numLagPts = 4;
- B.exposants = generateExposantsTetrahedron(o);
- B.points = gmshGeneratePointsTetrahedron(o,false);
- B.matrixLag2Bez = generateLag2BezMatrix(B.exposants, B.points, o, 3);
- std::vector< fullMatrix<double> > subPoints = generateSubPointsTetrahedron(o);
- B.divisor = generateDivisor(B.exposants, subPoints, B.matrixLag2Bez, o, 3);
- const polynomialBasis *F = polynomialBases::find(tag);
- generateGradShapes(B, B.points, F->monomials, F->coefficients);
- goto end;
- }
-
-
- switch (tag) {
- case MSH_QUA_4 : order = 1; break;
- case MSH_QUA_8 :
- case MSH_QUA_9 : order = 2; break;
- case MSH_QUA_12 :
- case MSH_QUA_16 : order = 3; break;
- case MSH_QUA_16I :
- case MSH_QUA_25 : order = 4; break;
- case MSH_QUA_20 :
- case MSH_QUA_36 : order = 5; break;
- case MSH_QUA_24 :
- case MSH_QUA_49 : order = 6; break;
- case MSH_QUA_28 :
- case MSH_QUA_64 : order = 7; break;
- case MSH_QUA_32 :
- case MSH_QUA_81 : order = 8; break;
- case MSH_QUA_36I :
- case MSH_QUA_100 : order = 9; break;
- case MSH_QUA_40 :
- case MSH_QUA_121 : order = 10; break;
- default: order = -1; break;
- }
- if (order >= 0) {
- int o = 2 * order - 1;
- B.numLagPts = 4;
- B.exposants = generateExposantsQuad(o);
- B.points = generatePointsQuad(o,false);
- B.matrixLag2Bez = generateLag2BezMatrix(B.exposants, B.points, o, 0);
- std::vector< fullMatrix<double> > subPoints = generateSubPointsQuad(o);
- B.divisor = generateDivisor(B.exposants, subPoints, B.matrixLag2Bez, o, 0);
- const polynomialBasis *F = polynomialBases::find(tag);
- generateGradShapes(B, B.points, F->monomials, F->coefficients);
- goto end;
- }
-
-
- switch (tag) {
- case MSH_PRI_6 : order = 1; break;
- case MSH_PRI_18 : order = 2; break;
- default: order = -1; break;
- }
- if (order >= 0) {
- int o = order * 3 - 1; // o=3*ord-2 on triangle base and =3*ord-1 for third dimension
- B.numLagPts = 4;
- B.exposants = generateExposantsPrism(o);
- B.points = generatePointsPrism(o,false);
- B.matrixLag2Bez = generateLag2BezMatrix(B.exposants, B.points, o, 2);
- std::vector< fullMatrix<double> > subPoints = generateSubPointsPrism(o);
- B.divisor = generateDivisor(B.exposants, subPoints, B.matrixLag2Bez, o, 2);
- const polynomialBasis *F = polynomialBases::find(tag);
- generateGradShapes(B, B.points, F->monomials, F->coefficients);
- }
- else {
- Msg::Error("Unknown function space %d: reverting to TET_4", tag);
- B.numLagPts = 4;
- B.exposants = generateExposantsTetrahedron(0);
- B.points = gmshGeneratePointsTetrahedron(0,false);
- B.matrixLag2Bez = generateLag2BezMatrix(B.exposants, B.points, 0, 3);
- std::vector< fullMatrix<double> > subPoints = generateSubPointsTetrahedron(0);
- B.divisor = generateDivisor(B.exposants, subPoints, B.matrixLag2Bez, 0, 3);
- const polynomialBasis *F = polynomialBases::find(tag);
- generateGradShapes(B, B.points, F->monomials, F->coefficients);
- }
-
-
-end :
-
- B.numDivisions = (int) pow(2., (int) B.points.size2());
-
- fs.insert(std::make_pair(tag, B));
- return &fs[tag];
-}
+// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#include "GmshDefines.h"
+#include "GmshMessage.h"
+#include "JacobianBasis.h"
+#include <vector>
+#include "polynomialBasis.h"
+
+// Bezier Exposants
+static fullMatrix<double> generate1DExposants(int order)
+{
+ fullMatrix<double> exposants(order + 1, 1);
+ exposants(0, 0) = 0;
+ if (order > 0) {
+ exposants(1, 0) = order;
+ for (int i = 2; i < order + 1; i++) exposants(i, 0) = i - 1;
+ }
+
+ return exposants;
+}
+
+static fullMatrix<double> generateExposantsTriangle(int order)
+{
+ int nbPoints = (order + 1) * (order + 2) / 2;
+ fullMatrix<double> exposants(nbPoints, 2);
+
+ exposants(0, 0) = 0;
+ exposants(0, 1) = 0;
+
+ if (order > 0) {
+ exposants(1, 0) = order;
+ exposants(1, 1) = 0;
+ exposants(2, 0) = 0;
+ exposants(2, 1) = order;
+
+ if (order > 1) {
+ int index = 3, ksi = 0, eta = 0;
+
+ for (int i = 0; i < order - 1; i++, index++) {
+ ksi++;
+ exposants(index, 0) = ksi;
+ exposants(index, 1) = eta;
+ }
+
+ ksi = order;
+
+ for (int i = 0; i < order - 1; i++, index++) {
+ ksi--;
+ eta++;
+ exposants(index, 0) = ksi;
+ exposants(index, 1) = eta;
+ }
+
+ eta = order;
+ ksi = 0;
+
+ for (int i = 0; i < order - 1; i++, index++) {
+ eta--;
+ exposants(index, 0) = ksi;
+ exposants(index, 1) = eta;
+ }
+
+ if (order > 2) {
+ fullMatrix<double> inner = generateExposantsTriangle(order - 3);
+ inner.add(1);
+ exposants.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
+ }
+ }
+ }
+
+ return exposants;
+}
+static fullMatrix<double> generateExposantsQuad(int order)
+{
+ int nbPoints = (order+1)*(order+1);
+ fullMatrix<double> exposants(nbPoints, 2);
+
+ exposants(0, 0) = 0;
+ exposants(0, 1) = 0;
+
+ if (order > 0) {
+ exposants(1, 0) = order;
+ exposants(1, 1) = 0;
+ exposants(2, 0) = order;
+ exposants(2, 1) = order;
+ exposants(3, 0) = 0;
+ exposants(3, 1) = order;
+
+ if (order > 1) {
+ int index = 4;
+ const static int edges[4][2]={{0,1},{1,2},{2,3},{3,0}};
+ for (int iedge=0; iedge<4; iedge++) {
+ int p0 = edges[iedge][0];
+ int p1 = edges[iedge][1];
+ for (int i = 1; i < order; i++, index++) {
+ exposants(index, 0) = exposants(p0, 0) + i*(exposants(p1,0)-exposants(p0,0))/order;
+ exposants(index, 1) = exposants(p0, 1) + i*(exposants(p1,1)-exposants(p0,1))/order;
+ }
+ }
+ if (order > 2) {
+ fullMatrix<double> inner = generateExposantsQuad(order - 2);
+ inner.add(1);
+ exposants.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
+ }
+ }
+ }
+
+ return exposants;
+}
+static int nbdoftriangle(int order) { return (order + 1) * (order + 2) / 2; }
+
+static void nodepositionface0(int order, double *u, double *v, double *w)
+{ // uv surface - orientation v0-v2-v1
+ int ndofT = nbdoftriangle(order);
+ if (order == 0) { u[0] = 0.; v[0] = 0.; w[0] = 0.; return; }
+
+ u[0]= 0.; v[0]= 0.; w[0] = 0.;
+ u[1]= 0.; v[1]= order; w[1] = 0.;
+ u[2]= order; v[2]= 0.; w[2] = 0.;
+
+ // 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] = k + 1;
+ v[3 + order - 1 + k] = order - 1 - k ;
+ w[3 + order - 1 + k] = 0.;
+
+ u[3 + 2 * (order - 1) + k] = order - 1 - k;
+ v[3 + 2 * (order - 1) + k] = 0.;
+ 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);
+ }
+ }
+ for (int k = 0; k < ndofT; k++){
+ u[k] = u[k] / order;
+ v[k] = v[k] / order;
+ w[k] = w[k] / order;
+ }
+}
+
+static void nodepositionface1(int order, double *u, double *v, double *w)
+{ // uw surface - orientation v0-v1-v3
+ int ndofT = nbdoftriangle(order);
+ if (order == 0) { u[0] = 0.; v[0] = 0.; w[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;
+ }
+ 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;
+ }
+}
+
+static void nodepositionface2(int order, double *u, double *v, double *w)
+{ // vw surface - orientation v0-v3-v2
+ int ndofT = nbdoftriangle(order);
+ if (order == 0) { u[0] = 0.; v[0] = 0.; return; }
+
+ u[0]= 0.; v[0]= 0.; w[0] = 0.;
+ u[1]= 0.; v[1]= 0.; w[1] = order;
+ u[2]= 0.; v[2]= order; w[2] = 0.;
+ // edges
+ for (int k = 0; k < (order - 1); k++){
+
+ u[3 + k] = 0.;
+ v[3 + k] = 0.;
+ w[3 + k] = k + 1;
+
+ u[3 + order - 1 + k] = 0.;
+ v[3 + order - 1 + k] = k + 1;
+ w[3 + order - 1 + k] = order - 1 - k;
+
+ 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);
+ 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);
+ 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;
+ }
+}
+
+static void nodepositionface3(int order, double *u, double *v, double *w)
+{ // uvw surface - orientation v3-v1-v2
+ int ndofT = nbdoftriangle(order);
+ if (order == 0) { u[0] = 0.; v[0] = 0.; w[0] = 0.; return; }
+
+ u[0]= 0.; v[0]= 0.; w[0] = order;
+ 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] = order - 1 - k;
+
+ 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] = k + 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) + 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;
+ }
+}
+
+static fullMatrix<double> generateExposantsTetrahedron(int order)
+{
+ int nbPoints = (order + 1) * (order + 2) * (order + 3) / 6;
+
+ fullMatrix<double> exposants(nbPoints, 3);
+
+ exposants(0, 0) = 0;
+ exposants(0, 1) = 0;
+ exposants(0, 2) = 0;
+
+ if (order > 0) {
+ exposants(1, 0) = order;
+ exposants(1, 1) = 0;
+ exposants(1, 2) = 0;
+
+ exposants(2, 0) = 0;
+ exposants(2, 1) = order;
+ exposants(2, 2) = 0;
+
+ exposants(3, 0) = 0;
+ exposants(3, 1) = 0;
+ exposants(3, 2) = order;
+
+ // edges e5 and e6 switched in original version, opposite direction
+ // the template has been defined in table edges_tetra and faces_tetra (MElement.h)
+
+ if (order > 1) {
+ for (int k = 0; k < (order - 1); k++) {
+ exposants(4 + k, 0) = k + 1;
+ exposants(4 + order - 1 + k, 0) = order - 1 - k;
+ exposants(4 + 2 * (order - 1) + k, 0) = 0;
+ exposants(4 + 3 * (order - 1) + k, 0) = 0;
+ // exposants(4 + 4 * (order - 1) + k, 0) = order - 1 - k;
+ // exposants(4 + 5 * (order - 1) + k, 0) = 0.;
+ exposants(4 + 4 * (order - 1) + k, 0) = 0;
+ exposants(4 + 5 * (order - 1) + k, 0) = k+1;
+
+ exposants(4 + k, 1) = 0;
+ exposants(4 + order - 1 + k, 1) = k + 1;
+ exposants(4 + 2 * (order - 1) + k, 1) = order - 1 - k;
+ exposants(4 + 3 * (order - 1) + k, 1) = 0;
+ // exposants(4 + 4 * (order - 1) + k, 1) = 0.;
+ // exposants(4 + 5 * (order - 1) + k, 1) = order - 1 - k;
+ exposants(4 + 4 * (order - 1) + k, 1) = k+1;
+ exposants(4 + 5 * (order - 1) + k, 1) = 0;
+
+ exposants(4 + k, 2) = 0;
+ exposants(4 + order - 1 + k, 2) = 0;
+ exposants(4 + 2 * (order - 1) + k, 2) = 0;
+ exposants(4 + 3 * (order - 1) + k, 2) = order - 1 - k;
+ exposants(4 + 4 * (order - 1) + k, 2) = order - 1 - k;
+ exposants(4 + 5 * (order - 1) + k, 2) = order - 1 - k;
+ }
+
+ if (order > 2) {
+ 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);
+
+ // u-v plane
+
+ for (int i = 0; i < nbdofface; i++){
+ exposants(ns + i, 0) = u[i] * (order - 3) + 1;
+ exposants(ns + i, 1) = v[i] * (order - 3) + 1;
+ exposants(ns + i, 2) = w[i] * (order - 3);
+ }
+
+ ns = ns + nbdofface;
+
+ // u-w plane
+
+ nodepositionface1(order - 3, u, v, w);
+
+ for (int i=0; i < nbdofface; i++){
+ exposants(ns + i, 0) = u[i] * (order - 3) + 1;
+ exposants(ns + i, 1) = v[i] * (order - 3) ;
+ exposants(ns + i, 2) = w[i] * (order - 3) + 1;
+ }
+
+ // v-w plane
+
+ ns = ns + nbdofface;
+
+ nodepositionface2(order - 3, u, v, w);
+
+ for (int i = 0; i < nbdofface; i++){
+ exposants(ns + i, 0) = u[i] * (order - 3);
+ exposants(ns + i, 1) = v[i] * (order - 3) + 1;
+ exposants(ns + i, 2) = w[i] * (order - 3) + 1;
+ }
+
+ // u-v-w plane
+
+ ns = ns + nbdofface;
+
+ nodepositionface3(order - 3, u, v, w);
+
+ for (int i = 0; i < nbdofface; i++){
+ exposants(ns + i, 0) = u[i] * (order - 3) + 1;
+ exposants(ns + i, 1) = v[i] * (order - 3) + 1;
+ exposants(ns + i, 2) = w[i] * (order - 3) + 1;
+ }
+
+ ns = ns + nbdofface;
+
+ delete [] u;
+ delete [] v;
+ delete [] w;
+
+ if (order > 3) {
+
+ fullMatrix<double> interior = generateExposantsTetrahedron(order - 4);
+ for (int k = 0; k < interior.size1(); k++) {
+ exposants(ns + k, 0) = 1 + interior(k, 0);
+ exposants(ns + k, 1) = 1 + interior(k, 1);
+ exposants(ns + k, 2) = 1 + interior(k, 2);
+ }
+ }
+ }
+ }
+ }
+
+ return exposants;
+}
+
+static fullMatrix<double> generateExposantsPrism(int order)
+{
+// int nbMonomials = (order + 1) * (order + 1) * (order + 2) / 2;
+//
+// fullMatrix<double> monomials(nbMonomials, 3);
+// int index = 0;
+// fullMatrix<double> lineMonoms = generate1DMonomials(order);
+// fullMatrix<double> triMonoms = generateExposantsTriangle(order);
+// // store monomials in right order
+// for (int currentOrder = 0; currentOrder <= order; currentOrder++) {
+// int orderT = currentOrder, orderL = currentOrder;
+// for (orderL = 0; orderL < currentOrder; orderL++) {
+// // do all permutations of monoms for orderL, orderT
+// int iL = orderL;
+// for (int iT = (orderT)*(orderT+1)/2; iT < (orderT+1)*(orderT+2)/2 ;iT++) {
+// monomials(index,0) = triMonoms(iT,0);
+// monomials(index,1) = triMonoms(iT,1);
+// monomials(index,2) = lineMonoms(iL,0);
+// index ++;
+// }
+// }
+// orderL = currentOrder;
+// for (orderT = 0; orderT <= currentOrder; orderT++) {
+// int iL = orderL;
+// for (int iT = (orderT)*(orderT+1)/2; iT < (orderT+1)*(orderT+2)/2 ;iT++) {
+// monomials(index,0) = triMonoms(iT,0);
+// monomials(index,1) = triMonoms(iT,1);
+// monomials(index,2) = lineMonoms(iL,0);
+// index ++;
+// }
+// }
+// }
+//// monomials.print("Pri monoms");
+// return monomials;
+
+ //const double prism18Pts[18][3] = {
+ // {0, 0, -1}, // 0
+ // {1, 0, -1}, // 1
+ // {0, 1, -1}, // 2
+ // {0, 0, 1}, // 3
+ // {1, 0, 1}, // 4
+ // {0, 1, 1}, // 5
+ // {0.5, 0, -1}, // 6
+ // {0, 0.5, -1}, // 7
+ // {0, 0, 0}, // 8
+ // {0.5, 0.5, -1}, // 9
+ // {1, 0, 0}, // 10
+ // {0, 1, 0}, // 11
+ // {0.5, 0, 1}, // 12
+ // {0, 0.5, 1}, // 13
+ // {0.5, 0.5, 1}, // 14
+ // {0.5, 0, 0}, // 15
+ // {0, 0.5, 0}, // 16
+ // {0.5, 0.5, 0}, // 17
+ // };
+
+ int nbPoints = (order + 1)*(order + 1)*(order + 2)/2;
+ fullMatrix<double> exposants(nbPoints, 3);
+
+ int index = 0;
+ fullMatrix<double> triExp = generateExposantsTriangle(order);
+ fullMatrix<double> lineExp = generate1DExposants(order);
+
+ /* if (order == 2)
+ for (int i =0; i<18; i++)
+ for (int j=0; j<3;j++)
+ exposants(i,j) = prism18Pts[i][j];
+ else*/
+ {
+ for (int i = 0; i < 3; i++) {
+ exposants(index,0) = triExp(i,0);
+ exposants(index,1) = triExp(i,1);
+ exposants(index,2) = lineExp(0,0);
+ index ++;
+ exposants(index,0) = triExp(i,0);
+ exposants(index,1) = triExp(i,1);
+ exposants(index,2) = lineExp(1,0);
+ index ++;
+ }
+ for (int i = 3; i < triExp.size1(); i++) {
+ exposants(index,0) = triExp(i,0);
+ exposants(index,1) = triExp(i,1);
+ exposants(index,2) = lineExp(0,0);
+ index ++;
+ exposants(index,0) = triExp(i,0);
+ exposants(index,1) = triExp(i,1);
+ exposants(index,2) = lineExp(1,0);
+ index ++;
+ }
+ for (int j = 2; j <lineExp.size1() ; j++) {
+ for (int i = 0; i < triExp.size1(); i++) {
+ exposants(index,0) = triExp(i,0);
+ exposants(index,1) = triExp(i,1);
+ exposants(index,2) = lineExp(j,0);
+ index ++;
+ }
+ }
+ }
+
+ return exposants;
+}
+
+// Sampling Points
+static fullMatrix<double> generate1DPoints(int order)
+{
+ fullMatrix<double> line(order + 1, 1);
+ line(0,0) = 0;
+ if (order > 0) {
+ line(0, 0) = 0.;
+ line(1, 0) = 1.;
+ double dd = 1. / order;
+ for (int i = 2; i < order + 1; i++) line(i, 0) = dd * (i - 1);
+ }
+
+ return line;
+}
+
+static fullMatrix<double> 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);
+
+ fullMatrix<double> point(nbPoints, 3);
+
+ double overOrder = (order == 0 ? 1. : 1. / order);
+
+ 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(2, 0) = 0.;
+ point(2, 1) = order;
+ point(2, 2) = 0.;
+
+ point(3, 0) = 0.;
+ point(3, 1) = 0.;
+ point(3, 2) = order;
+
+ // edges e5 and e6 switched in original version, opposite direction
+ // the template has been defined in table edges_tetra and faces_tetra (MElement.h)
+
+ 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.;
+ point(4 + 4 * (order - 1) + k, 0) = 0.;
+ point(4 + 5 * (order - 1) + k, 0) = k+1;
+
+ 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 + 4 * (order - 1) + k, 1) = k+1;
+ point(4 + 5 * (order - 1) + k, 1) = 0.;
+
+ 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) = order - 1 - k;
+ point(4 + 4 * (order - 1) + k, 2) = order - 1 - k;
+ point(4 + 5 * (order - 1) + k, 2) = order - 1 - k;
+ }
+
+ if (order > 2) {
+ 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);
+
+ // u-v plane
+
+ 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;
+
+ // u-w plane
+
+ 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.;
+ }
+
+ // v-w plane
+
+ ns = ns + nbdofface;
+
+ nodepositionface2(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.;
+ }
+
+ // u-v-w plane
+
+ ns = ns + nbdofface;
+
+ nodepositionface3(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.;
+ }
+
+ ns = ns + nbdofface;
+
+ delete [] u;
+ delete [] v;
+ delete [] w;
+
+ if (!serendip && order > 3) {
+
+ fullMatrix<double> interior = gmshGeneratePointsTetrahedron(order - 4, false);
+ for (int 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);
+ }
+ }
+ }
+ }
+ }
+
+ point.scale(overOrder);
+ return point;
+}
+
+static fullMatrix<double> gmshGeneratePointsTriangle(int order, bool serendip)
+{
+ int nbPoints = serendip ? 3 * order : (order + 1) * (order + 2) / 2;
+ fullMatrix<double> point(nbPoints, 2);
+
+ point(0, 0) = 0;
+ point(0, 1) = 0;
+
+ if (order > 0) {
+ double dd = 1. / order;
+
+ point(1, 0) = 1;
+ point(1, 1) = 0;
+ point(2, 0) = 0;
+ point(2, 1) = 1;
+
+ int index = 3;
+
+ if (order > 1) {
+
+ double ksi = 0;
+ double eta = 0;
+
+ for (int i = 0; i < order - 1; i++, index++) {
+ ksi += dd;
+ point(index, 0) = ksi;
+ point(index, 1) = eta;
+ }
+
+ ksi = 1.;
+
+ for (int i = 0; i < order - 1; i++, index++) {
+ ksi -= dd;
+ eta += dd;
+ point(index, 0) = ksi;
+ point(index, 1) = eta;
+ }
+
+ eta = 1.;
+ ksi = 0.;
+
+ for (int i = 0; i < order - 1; i++, index++) {
+ eta -= dd;
+ point(index, 0) = ksi;
+ point(index, 1) = eta;
+ }
+
+ if (order > 2 && !serendip) {
+ fullMatrix<double> inner = gmshGeneratePointsTriangle(order - 3, serendip);
+ inner.scale(1. - 3. * dd);
+ inner.add(dd);
+ point.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
+ }
+ }
+ }
+ return point;
+}
+
+static fullMatrix<double> generatePointsPrism(int order, bool serendip)
+{
+ const double prism18Pts[18][3] = {
+ {0, 0, 0}, // 0
+ {1, 0, 0}, // 1
+ {0, 1, 0}, // 2
+ {0, 0, 1}, // 3
+ {1, 0, 1}, // 4
+ {0, 1, 1}, // 5
+ {.5, 0, 0}, // 6
+ {0, .5, 0}, // 7
+ {0, 0, .5}, // 8
+ {.5, .5, 0}, // 9
+ {1, 0, .5}, // 10
+ {0, 1, .5}, // 11
+ {.5, 0, 1}, // 12
+ {0, .5, 1}, // 13
+ {.5, .5, 1}, // 14
+ {.5, 0, .5}, // 15
+ {0, .5, .5}, // 16
+ {.5, .5, .5}, // 17
+ };
+
+ int nbPoints = (order + 1)*(order + 1)*(order + 2)/2;
+ fullMatrix<double> point(nbPoints, 3);
+
+ int index = 0;
+
+ if (order == 2)
+ for (int i =0; i<18; i++)
+ for (int j=0; j<3;j++)
+ point(i,j) = prism18Pts[i][j];
+ else {
+ fullMatrix<double> triPoints = gmshGeneratePointsTriangle(order,false);
+ fullMatrix<double> linePoints = generate1DPoints(order);
+ for (int j = 0; j <linePoints.size1() ; j++) {
+ for (int i = 0; i < triPoints.size1(); i++) {
+ point(index,0) = triPoints(i,0);
+ point(index,1) = triPoints(i,1);
+ point(index,2) = linePoints(j,0);
+ index ++;
+ }
+ }
+ }
+ return point;
+}
+
+static fullMatrix<double> generatePointsQuad(int order, bool serendip)
+{
+ int nbPoints = serendip ? order*4 : (order+1)*(order+1);
+ fullMatrix<double> point(nbPoints, 2);
+
+ if (order > 0) {
+ point(0, 0) = 0;
+ point(0, 1) = 0;
+ point(1, 0) = 1;
+ point(1, 1) = 0;
+ point(2, 0) = 1;
+ point(2, 1) = 1;
+ point(3, 0) = 0;
+ point(3, 1) = 1;
+
+ if (order > 1) {
+ int index = 4;
+ const static int edges[4][2]={{0,1},{1,2},{2,3},{3,0}};
+ for (int iedge=0; iedge<4; iedge++) {
+ int p0 = edges[iedge][0];
+ int p1 = edges[iedge][1];
+ for (int i = 1; i < order; i++, index++) {
+ point(index, 0) = point(p0, 0) + i*(point(p1,0)-point(p0,0))/order;
+ point(index, 1) = point(p0, 1) + i*(point(p1,1)-point(p0,1))/order;
+ }
+ }
+ if (order > 2 && !serendip) {
+ fullMatrix<double> inner = generatePointsQuad(order - 2, false);
+ inner.scale(1. - 2./order);
+ point.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
+ }
+ }
+ }
+ else {
+ point(0, 0) = 0;
+ point(0, 1) = 0;
+ }
+ return point;
+}
+
+// Sub Control Points
+static std::vector< fullMatrix<double> > generateSubPointsLine(int order)
+{
+ std::vector< fullMatrix<double> > subPoints(2);
+ fullMatrix<double> prox;
+
+ subPoints[0] = generate1DPoints(order);
+ subPoints[0].scale(.5);
+
+ subPoints[1].copy(subPoints[0]);
+ subPoints[1].add(.5);
+
+ return subPoints;
+}
+
+static std::vector< fullMatrix<double> > generateSubPointsTriangle(int order)
+{
+ std::vector< fullMatrix<double> > subPoints(4);
+ fullMatrix<double> prox;
+
+ subPoints[0] = gmshGeneratePointsTriangle(order, false);
+ subPoints[0].scale(.5);
+
+ subPoints[1].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[1], 0, 1);
+ prox.add(.5);
+
+ subPoints[2].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[2], 1, 1);
+ prox.add(.5);
+
+ subPoints[3].copy(subPoints[0]);
+ subPoints[3].scale(-1.);
+ subPoints[3].add(.5);
+
+ return subPoints;
+}
+
+static std::vector< fullMatrix<double> > generateSubPointsTetrahedron(int order)
+{
+ std::vector< fullMatrix<double> > subPoints(8);
+ fullMatrix<double> prox1;
+ fullMatrix<double> prox2;
+
+ subPoints[0] = gmshGeneratePointsTetrahedron(order, false);
+ subPoints[0].scale(.5);
+
+ subPoints[1].copy(subPoints[0]);
+ prox1.setAsProxy(subPoints[1], 0, 1);
+ prox1.add(.5);
+
+ subPoints[2].copy(subPoints[0]);
+ prox1.setAsProxy(subPoints[2], 1, 1);
+ prox1.add(.5);
+
+ subPoints[3].copy(subPoints[0]);
+ prox1.setAsProxy(subPoints[3], 2, 1);
+ prox1.add(.5);
+
+ // u := .5-u-w
+ // v := .5-v-w
+ // w := w
+ subPoints[4].copy(subPoints[0]);
+ prox1.setAsProxy(subPoints[4], 0, 2);
+ prox1.scale(-1.);
+ prox1.add(.5);
+ prox1.setAsProxy(subPoints[4], 0, 1);
+ prox2.setAsProxy(subPoints[4], 2, 1);
+ prox1.add(prox2, -1.);
+ prox1.setAsProxy(subPoints[4], 1, 1);
+ prox1.add(prox2, -1.);
+
+ // u := u
+ // v := .5-v
+ // w := w+v
+ subPoints[5].copy(subPoints[0]);
+ prox1.setAsProxy(subPoints[5], 2, 1);
+ prox2.setAsProxy(subPoints[5], 1, 1);
+ prox1.add(prox2);
+ prox2.scale(-1.);
+ prox2.add(.5);
+
+ // u := .5-u
+ // v := v
+ // w := w+u
+ subPoints[6].copy(subPoints[0]);
+ prox1.setAsProxy(subPoints[6], 2, 1);
+ prox2.setAsProxy(subPoints[6], 0, 1);
+ prox1.add(prox2);
+ prox2.scale(-1.);
+ prox2.add(.5);
+
+ // u := u+w
+ // v := v+w
+ // w := .5-w
+ subPoints[7].copy(subPoints[0]);
+ prox1.setAsProxy(subPoints[7], 0, 1);
+ prox2.setAsProxy(subPoints[7], 2, 1);
+ prox1.add(prox2);
+ prox1.setAsProxy(subPoints[7], 1, 1);
+ prox1.add(prox2);
+ prox2.scale(-1.);
+ prox2.add(.5);
+
+
+ return subPoints;
+}
+
+static std::vector< fullMatrix<double> > generateSubPointsQuad(int order)
+{
+ std::vector< fullMatrix<double> > subPoints(4);
+ fullMatrix<double> prox;
+
+ subPoints[0] = generatePointsQuad(order, false);
+ subPoints[0].scale(.5);
+
+ subPoints[1].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[1], 0, 1);
+ prox.add(.5);
+
+ subPoints[2].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[2], 1, 1);
+ prox.add(.5);
+
+ subPoints[3].copy(subPoints[1]);
+ prox.setAsProxy(subPoints[3], 1, 1);
+ prox.add(.5);
+
+ return subPoints;
+}
+
+static std::vector< fullMatrix<double> > generateSubPointsPrism(int order)
+{
+ std::vector< fullMatrix<double> > subPoints(8);
+ fullMatrix<double> prox;
+
+ subPoints[0] = generatePointsPrism(order, false);
+ subPoints[0].scale(.5);
+
+ subPoints[1].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[1], 0, 1);
+ prox.add(.5);
+
+ subPoints[2].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[2], 1, 1);
+ prox.add(.5);
+
+ subPoints[3].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[3], 0, 2);
+ prox.scale(-1.);
+ prox.add(.5);
+
+ subPoints[4].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[4], 2, 1);
+ prox.add(.5);
+
+ subPoints[5].copy(subPoints[1]);
+ prox.setAsProxy(subPoints[5], 2, 1);
+ prox.add(.5);
+
+ subPoints[6].copy(subPoints[2]);
+ prox.setAsProxy(subPoints[6], 2, 1);
+ prox.add(.5);
+
+ subPoints[7].copy(subPoints[3]);
+ prox.setAsProxy(subPoints[7], 2, 1);
+ prox.add(.5);
+
+ return subPoints;
+}
+
+// Matrices generation
+static int nChoosek(int n, int k)
+{
+ if (n < k || k < 0) {
+ Msg::Error("Wrong argument for combination.");
+ return 1;
+ }
+
+ if (k > n/2) k = n-k;
+ if (k == 1)
+ return n;
+ if (k == 0)
+ return 1;
+
+ int c = 1;
+ for (int i = 1; i <= k; i++, n--) (c *= n) /= i;
+ return c;
+}
+
+
+static fullMatrix<double> generateLag2BezMatrix
+ (const fullMatrix<double> &exposant, const fullMatrix<double> &point,
+ int order, int dimSimplex, bool invert = true)
+{
+
+ if(exposant.size1() != point.size1() || exposant.size2() != point.size2()){
+ Msg::Fatal("Wrong sizes for Bezier coefficients generation %d %d -- %d %d",
+ exposant.size1(),point.size1(),
+ exposant.size2(),point.size2());
+ return fullMatrix<double>(1, 1);
+ }
+
+ int ndofs = exposant.size1();
+ int dim = exposant.size2();
+
+ fullMatrix<double> Vandermonde(ndofs, ndofs);
+ for (int i = 0; i < ndofs; i++) {
+ for (int j = 0; j < ndofs; j++) {
+ double dd = 1.;
+
+ double pointCompl = 1.;
+ int exposantCompl = order;
+ for (int k = 0; k < dimSimplex; k++) {
+ dd *= nChoosek(exposantCompl, (int) exposant(i, k))
+ * pow(point(j, k), exposant(i, k));
+ pointCompl -= point(j, k);
+ exposantCompl -= (int) exposant(i, k);
+ }
+ dd *= pow(pointCompl, exposantCompl);
+
+ for (int k = dimSimplex; k < dim; k++)
+ dd *= nChoosek(order, (int) exposant(i, k))
+ * pow(point(j, k), exposant(i, k))
+ * pow(1. - point(j, k), order - exposant(i, k));
+
+ Vandermonde(j, i) = dd;
+ }
+ }
+
+ if (!invert) return Vandermonde;
+
+ fullMatrix<double> coefficient(ndofs, ndofs);
+ Vandermonde.invert(coefficient);
+ return coefficient;
+}
+
+
+
+static fullMatrix<double> generateDivisor
+ (const fullMatrix<double> &exposants, const std::vector< fullMatrix<double> > &subPoints,
+ const fullMatrix<double> &lag2Bez, int order, int dimSimplex)
+{
+ if (exposants.size1() != lag2Bez.size1() || exposants.size1() != lag2Bez.size2()){
+ Msg::Fatal("Wrong sizes for Bezier Divisor %d %d -- %d %d",
+ exposants.size1(), lag2Bez.size1(),
+ exposants.size1(), lag2Bez.size2());
+ return fullMatrix<double>(1, 1);
+ }
+
+ int nbPts = lag2Bez.size1();
+ int nbSubPts = nbPts * subPoints.size();
+
+ fullMatrix<double> intermediate2(nbPts, nbPts);
+ fullMatrix<double> divisor(nbSubPts, nbPts);
+
+ for (unsigned int i = 0; i < subPoints.size(); i++) {
+ fullMatrix<double> intermediate1 =
+ generateLag2BezMatrix(exposants, subPoints[i], order, dimSimplex, false);
+ lag2Bez.mult(intermediate1, intermediate2);
+ divisor.copy(intermediate2, 0, nbPts, 0, nbPts, i*nbPts, 0);
+ }
+ return divisor;
+}
+
+
+
+static void generateGradShapes(JacobianBasis &jfs, const fullMatrix<double> &points
+ , const fullMatrix<double> &monomials, const fullMatrix<double> &coefficients)
+{
+
+ /*{
+ Msg::Info("Printing vector jacobian':");
+ int ni = points.size1();
+ int nj = points.size2();
+ Msg::Info(" ");
+ for(int I = 0; I < ni; I++){
+ Msg::Info("%lf - %lf", points(I, 0), points(I, 1));
+ }
+ Msg::Info(" ");
+ }
+ {
+ Msg::Info("Printing vector jacobian':");
+ int ni = monomials.size1();
+ int nj = monomials.size2();
+ Msg::Info(" ");
+ for(int I = 0; I < ni; I++){
+ Msg::Info("%lf - %lf", monomials(I, 0), monomials(I, 1));
+ }
+ Msg::Info(" ");
+ }
+ {
+ Msg::Info("Printing vector jacobian':");
+ int ni = coefficients.size1();
+ int nj = coefficients.size2();
+ Msg::Info(" ");
+ for(int I = 0; I < ni; I++){
+ Msg::Info(" ");
+ for(int J = 0; J < nj; J++){
+ Msg::Info("%lf", coefficients(J, I));
+ }
+ }
+ Msg::Info(" ");
+ }*/
+
+ int nbPts = points.size1();
+ int nbDofs = monomials.size1();
+ int dim = points.size2();
+
+ switch (dim) {
+ case 3 :
+ jfs.gradShapeMatZ.resize(nbPts, nbDofs, true);
+ case 2 :
+ jfs.gradShapeMatY.resize(nbPts, nbDofs, true);
+ case 1 :
+ jfs.gradShapeMatX.resize(nbPts, nbDofs, true);
+ break;
+ default :
+ return;
+ }
+
+ double dx, dy, dz;
+
+ switch (dim) {
+ case 3 :
+ for (int i = 0; i < nbDofs; i++) {
+ for (int k = 0; k < nbPts; k++) {
+
+ if ((int) monomials(i, 0) > 0) {
+ dx = pow( points(k, 0), monomials(i, 0)-1 ) * monomials(i, 0)
+ * pow( points(k, 1), monomials(i, 1) )
+ * pow( points(k, 2), monomials(i, 2) );
+ for (int j = 0; j < nbDofs; j++)
+ jfs.gradShapeMatX(k, j) += coefficients(j, i) * dx;
+ }
+ if ((int) monomials(i, 1) > 0.) {
+ dy = pow( points(k, 0), monomials(i, 0) )
+ * pow( points(k, 1), monomials(i, 1)-1 ) * monomials(i, 1)
+ * pow( points(k, 2), monomials(i, 2) );
+ for (int j = 0; j < nbDofs; j++)
+ jfs.gradShapeMatY(k, j) += coefficients(j, i) * dy;
+ }
+ if ((int) monomials(i, 2) > 0.) {
+ dz = pow( points(k, 0), monomials(i, 0) )
+ * pow( points(k, 1), monomials(i, 1) )
+ * pow( points(k, 2), monomials(i, 2)-1 ) * monomials(i, 2);
+ for (int j = 0; j < nbDofs; j++)
+ jfs.gradShapeMatZ(k, j) += coefficients(j, i) * dz;
+ }
+ }
+ }
+ return;
+
+ case 2 :
+ for (int i = 0; i < nbDofs; i++) {
+ for (int k = 0; k < nbPts; k++) {
+
+ if ((int) monomials(i, 0) > 0) {
+ dx = pow( points(k, 0), (int) monomials(i, 0)-1 ) * monomials(i, 0)
+ * pow( points(k, 1), (int) monomials(i, 1) );
+ for (int j = 0; j < nbDofs; j++)
+ jfs.gradShapeMatX(k, j) += coefficients(j, i) * dx;
+ }
+ if ((int) monomials(i, 1) > 0) {
+ dy = pow( points(k, 0), (int) monomials(i, 0) )
+ * pow( points(k, 1), (int) monomials(i, 1)-1 ) * monomials(i, 1);
+ for (int j = 0; j < nbDofs; j++)
+ jfs.gradShapeMatY(k, j) += coefficients(j, i) * dy;
+ }
+ }
+ }
+ return;
+
+ case 1 :
+ for (int i = 0; i < nbDofs; i++) {
+ for (int k = 0; k < nbPts; k++) {
+
+ if ((int) monomials(i, 0) > 0) {
+ dx = pow( points(k, 0), (int) monomials(i, 0)-1 ) * monomials(i, 0);
+ for (int j = 0; j < nbDofs; j++)
+ jfs.gradShapeMatX(k, j) += coefficients(j, i) * dx;
+ }
+ }
+ }
+ }
+ return;
+}
+
+std::map<int, JacobianBasis> JacobianBases::fs;
+
+const JacobianBasis *JacobianBases::find(int tag)
+{
+ std::map<int, JacobianBasis>::const_iterator it = fs.find(tag);
+ if (it != fs.end()) return &it->second;
+
+ JacobianBasis B;
+ B.numLagPts = -1;
+
+ int order;
+
+
+ if (tag == MSH_PNT) {
+ B.numLagPts = 1;
+ B.exposants = generate1DExposants(0);
+ B.points = generate1DPoints(0);
+ B.matrixLag2Bez = generateLag2BezMatrix(B.exposants, B.points, 0, 0);
+ /*std::vector< fullMatrix<double> > subPoints = generateSubPointsLine(0);
+ B.divisor = generateDivisor(B.exposants, subPoints, B.matrixLag2Bez, 0, 0);
+ const polynomialBasis *F = polynomialBases::find(tag);
+ generateGradShapes(B, B.points, F->monomials, F->coefficients);*/
+ goto end;
+ }
+
+
+ switch (tag) {
+ case MSH_LIN_2: order = 1; break;
+ case MSH_LIN_3: order = 2; break;
+ case MSH_LIN_4: order = 3; break;
+ case MSH_LIN_5: order = 4; break;
+ case MSH_LIN_6: order = 5; break;
+ case MSH_LIN_7: order = 6; break;
+ case MSH_LIN_8: order = 7; break;
+ case MSH_LIN_9: order = 8; break;
+ case MSH_LIN_10: order = 9; break;
+ case MSH_LIN_11: order = 10; break;
+ default: order = -1; break;
+ }
+ if (order >= 0) {
+ int o = order - 1;
+ B.numLagPts = 2;
+ B.exposants = generate1DExposants(o);
+ B.points = generate1DPoints(o);
+ B.matrixLag2Bez = generateLag2BezMatrix(B.exposants, B.points, o, 0);
+ std::vector< fullMatrix<double> > subPoints = generateSubPointsLine(0);
+ B.divisor = generateDivisor(B.exposants, subPoints, B.matrixLag2Bez, o, 0);
+ const polynomialBasis *F = polynomialBases::find(tag);
+ generateGradShapes(B, B.points, F->monomials, F->coefficients);
+ goto end;
+ }
+
+
+
+ switch (tag) {
+ case MSH_TRI_3 : order = 1; break;
+ case MSH_TRI_6 : order = 2; break;
+ case MSH_TRI_9 :
+ case MSH_TRI_10 : order = 3; break;
+ case MSH_TRI_12 :
+ case MSH_TRI_15 : order = 4; break;
+ case MSH_TRI_15I :
+ case MSH_TRI_21 : order = 5; break;
+ case MSH_TRI_18 :
+ case MSH_TRI_28 : order = 6; break;
+ case MSH_TRI_21I :
+ case MSH_TRI_36 : order = 7; break;
+ case MSH_TRI_24 :
+ case MSH_TRI_45 : order = 8; break;
+ case MSH_TRI_27 :
+ case MSH_TRI_55 : order = 9; break;
+ case MSH_TRI_30 :
+ case MSH_TRI_66 : order = 10; break;
+ default: order = -1; break;
+ }
+ if (order >= 0) {
+ int o = 2 * (order - 1);
+ B.numLagPts = 3;
+ B.exposants = generateExposantsTriangle(o);
+ B.points = gmshGeneratePointsTriangle(o,false);
+ B.matrixLag2Bez = generateLag2BezMatrix(B.exposants, B.points, o, 2);
+ std::vector< fullMatrix<double> > subPoints = generateSubPointsTriangle(o);
+ B.divisor = generateDivisor(B.exposants, subPoints, B.matrixLag2Bez, o, 2);
+ const polynomialBasis *F = polynomialBases::find(tag);
+ generateGradShapes(B, B.points, F->monomials, F->coefficients);
+ goto end;
+ }
+
+
+ switch (tag) {
+ case MSH_TET_4 : order = 1; break;
+ case MSH_TET_10 : order = 2; break;
+ case MSH_TET_20 : order = 3; break;
+ case MSH_TET_34 :
+ case MSH_TET_35 : order = 4; break;
+ case MSH_TET_52 :
+ case MSH_TET_56 : order = 5; break;
+ case MSH_TET_84 : order = 6; break;
+ case MSH_TET_120 : order = 7; break;
+ case MSH_TET_165 : order = 8; break;
+ case MSH_TET_220 : order = 9; break;
+ case MSH_TET_286 : order = 10; break;
+ default: order = -1; break;
+ }
+ if (order >= 0) {
+ int o = 3 * (order - 1);
+ B.numLagPts = 4;
+ B.exposants = generateExposantsTetrahedron(o);
+ B.points = gmshGeneratePointsTetrahedron(o,false);
+ B.matrixLag2Bez = generateLag2BezMatrix(B.exposants, B.points, o, 3);
+ std::vector< fullMatrix<double> > subPoints = generateSubPointsTetrahedron(o);
+ B.divisor = generateDivisor(B.exposants, subPoints, B.matrixLag2Bez, o, 3);
+ const polynomialBasis *F = polynomialBases::find(tag);
+ generateGradShapes(B, B.points, F->monomials, F->coefficients);
+ goto end;
+ }
+
+
+ switch (tag) {
+ case MSH_QUA_4 : order = 1; break;
+ case MSH_QUA_8 :
+ case MSH_QUA_9 : order = 2; break;
+ case MSH_QUA_12 :
+ case MSH_QUA_16 : order = 3; break;
+ case MSH_QUA_16I :
+ case MSH_QUA_25 : order = 4; break;
+ case MSH_QUA_20 :
+ case MSH_QUA_36 : order = 5; break;
+ case MSH_QUA_24 :
+ case MSH_QUA_49 : order = 6; break;
+ case MSH_QUA_28 :
+ case MSH_QUA_64 : order = 7; break;
+ case MSH_QUA_32 :
+ case MSH_QUA_81 : order = 8; break;
+ case MSH_QUA_36I :
+ case MSH_QUA_100 : order = 9; break;
+ case MSH_QUA_40 :
+ case MSH_QUA_121 : order = 10; break;
+ default: order = -1; break;
+ }
+ if (order >= 0) {
+ int o = 2 * order - 1;
+ B.numLagPts = 4;
+ B.exposants = generateExposantsQuad(o);
+ B.points = generatePointsQuad(o,false);
+ B.matrixLag2Bez = generateLag2BezMatrix(B.exposants, B.points, o, 0);
+ std::vector< fullMatrix<double> > subPoints = generateSubPointsQuad(o);
+ B.divisor = generateDivisor(B.exposants, subPoints, B.matrixLag2Bez, o, 0);
+ const polynomialBasis *F = polynomialBases::find(tag);
+ generateGradShapes(B, B.points, F->monomials, F->coefficients);
+ goto end;
+ }
+
+
+ switch (tag) {
+ case MSH_PRI_6 : order = 1; break;
+ case MSH_PRI_18 : order = 2; break;
+ default: order = -1; break;
+ }
+ if (order >= 0) {
+ int o = order * 3 - 1; // o=3*ord-2 on triangle base and =3*ord-1 for third dimension
+ B.numLagPts = 4;
+ B.exposants = generateExposantsPrism(o);
+ B.points = generatePointsPrism(o,false);
+ B.matrixLag2Bez = generateLag2BezMatrix(B.exposants, B.points, o, 2);
+ std::vector< fullMatrix<double> > subPoints = generateSubPointsPrism(o);
+ B.divisor = generateDivisor(B.exposants, subPoints, B.matrixLag2Bez, o, 2);
+ const polynomialBasis *F = polynomialBases::find(tag);
+ generateGradShapes(B, B.points, F->monomials, F->coefficients);
+ }
+ else {
+ Msg::Error("Unknown function space %d: reverting to TET_4", tag);
+ B.numLagPts = 4;
+ B.exposants = generateExposantsTetrahedron(0);
+ B.points = gmshGeneratePointsTetrahedron(0,false);
+ B.matrixLag2Bez = generateLag2BezMatrix(B.exposants, B.points, 0, 3);
+ std::vector< fullMatrix<double> > subPoints = generateSubPointsTetrahedron(0);
+ B.divisor = generateDivisor(B.exposants, subPoints, B.matrixLag2Bez, 0, 3);
+ const polynomialBasis *F = polynomialBases::find(tag);
+ generateGradShapes(B, B.points, F->monomials, F->coefficients);
+ }
+
+
+end :
+
+ B.numDivisions = (int) pow(2., (int) B.points.size2());
+
+ fs.insert(std::make_pair(tag, B));
+ return &fs[tag];
+}
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index 645c512c7583567b54383b4dfeee9b740360be02..c3daaa23ecb35774f243690083e7100fe3e8bf34 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -1,35 +1,35 @@
-// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to <gmsh@geuz.org>.
-
-#ifndef _JACOBIAN_BASIS_H_
-#define _JACOBIAN_BASIS_H_
-
-#include <map>
-#include <vector>
-#include "fullMatrix.h"
-
-class JacobianBasis
-{
- public:
- int numLagPts;
- int numDivisions;
- fullMatrix<double> exposants; //exposants of Bezier FF
- fullMatrix<double> points; //sampling point
- fullMatrix<double> matrixLag2Bez;
- fullMatrix<double> gradShapeMatX;
- fullMatrix<double> gradShapeMatY;
- fullMatrix<double> gradShapeMatZ;
- fullMatrix<double> divisor;
-};
-
-class JacobianBases
-{
- private:
- static std::map<int, JacobianBasis> fs;
- public:
- static const JacobianBasis *find(int);
-};
-
-#endif
+// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _JACOBIAN_BASIS_H_
+#define _JACOBIAN_BASIS_H_
+
+#include <map>
+#include <vector>
+#include "fullMatrix.h"
+
+class JacobianBasis
+{
+ public:
+ int numLagPts;
+ int numDivisions;
+ fullMatrix<double> exposants; //exposants of Bezier FF
+ fullMatrix<double> points; //sampling point
+ fullMatrix<double> matrixLag2Bez;
+ fullMatrix<double> gradShapeMatX;
+ fullMatrix<double> gradShapeMatY;
+ fullMatrix<double> gradShapeMatZ;
+ fullMatrix<double> divisor;
+};
+
+class JacobianBases
+{
+ private:
+ static std::map<int, JacobianBasis> fs;
+ public:
+ static const JacobianBasis *find(int);
+};
+
+#endif
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 5659d11090cf029bd0af9d310fc241dac25a0ea3..8eb925445a2850959c6c66534d97f976179e4541 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -17,21 +17,22 @@
static void printClosure(polynomialBasis::clCont &fullClosure, std::vector<int> &closureRef,
polynomialBasis::clCont &closures)
{
- for (int i = 0; i <closures.size(); i++) {
+ for (unsigned int i = 0; i < closures.size(); i++) {
printf("%3i (%2i): ", i, closureRef[i]);
if(closureRef[i]==-1){
printf("\n");
continue;
}
- for (int j = 0; j < fullClosure[i].size(); j++) {
+ for (unsigned int j = 0; j < fullClosure[i].size(); j++) {
printf("%2i ", fullClosure[i][j]);
}
printf (" -- ");
- for (int j = 0; j < closures[closureRef[i]].size(); j++) {
+ for (unsigned int j = 0; j < closures[closureRef[i]].size(); j++) {
std::string equalSign = "-";
if (fullClosure[i][closures[closureRef[i]][j]] != closures[i][j])
equalSign = "#";
- printf("%2i%s%2i ", fullClosure[i][closures[closureRef[i]][j]],equalSign.c_str(), closures[i][j]);
+ printf("%2i%s%2i ", fullClosure[i][closures[closureRef[i]][j]],equalSign.c_str(),
+ closures[i][j]);
}
printf("\n");
}
@@ -51,7 +52,7 @@ static int getTriangleType (int order)
case 8 : return MSH_TRI_45;
case 9 : return MSH_TRI_55;
case 10 : return MSH_TRI_66;
- default : Msg::Error("triangle order %i unknown", order);
+ default : Msg::Error("triangle order %i unknown", order); return 0;
}
}
@@ -69,7 +70,7 @@ static int getQuadType (int order)
case 8 : return MSH_QUA_81;
case 9 : return MSH_QUA_100;
case 10 : return MSH_QUA_121;
- default : Msg::Error("quad order %i unknown", order);
+ default : Msg::Error("quad order %i unknown", order); return 0;
}
}
@@ -87,7 +88,7 @@ static int getLineType (int order)
case 8 : return MSH_LIN_9;
case 9 : return MSH_LIN_10;
case 10 : return MSH_LIN_11;
- default : Msg::Error("line order %i unknown", order);
+ default : Msg::Error("line order %i unknown", order); return 0;
}
}
@@ -162,7 +163,7 @@ const fullMatrix<double> &polynomialBasis::getGradientAtFaceIntegrationPoints(in
if (dfAtFace.empty()) {
dfAtFace.resize(closures.size());
int nbPsi = points.size1();
- for (int iClosure=0; iClosure< closures.size(); iClosure++) {
+ for (unsigned int iClosure = 0; iClosure < closures.size(); iClosure++) {
fullMatrix<double> integrationFace, weight;
const polynomialBasis *fsFace = polynomialBases::find(closures[iClosure].type);
gaussIntegration::get(fsFace->parentType, integrationOrder, integrationFace, weight);
@@ -170,7 +171,7 @@ const fullMatrix<double> &polynomialBasis::getGradientAtFaceIntegrationPoints(in
double f[1256];
for (int i = 0; i < integrationFace.size1(); i++){
fsFace->f(integrationFace(i,0),integrationFace(i,1),integrationFace(i,2),f);
- for(size_t j=0; j<closures[iClosure].size(); j++) {
+ for(unsigned int j = 0; j < closures[iClosure].size(); j++) {
int jNod = closures[iClosure][j];
for(int k = 0; k < points.size2(); k++)
integration(i,k) += f[j] * points(jNod,k);
@@ -179,8 +180,8 @@ const fullMatrix<double> &polynomialBasis::getGradientAtFaceIntegrationPoints(in
fullMatrix<double> &v = dfAtFace[iClosure];
v.resize(nbPsi, integration.size1()*3);
double g[1256][3];
- for (size_t xi = 0; xi< integration.size1(); xi++) {
- df(integration(xi,0 ), integration(xi,1), integration(xi,2), g);
+ for (int xi = 0; xi < integration.size1(); xi++) {
+ df(integration(xi, 0), integration(xi, 1), integration(xi, 2), g);
for (int j = 0; j < nbPsi; j++) {
v(j, xi*3) = g[j][0];
v(j, xi*3+1) = g[j][1];
@@ -1153,7 +1154,7 @@ static void addEdgeNodes(polynomialBasis::clCont &closureFull, const int *edges,
nodes2edges[edges[i]][edges[i + 1]] = i;
nodes2edges[edges[i + 1]][edges[i]] = i + 1;
}
- for (int iClosure = 0; iClosure < closureFull.size(); iClosure++) {
+ for (unsigned int iClosure = 0; iClosure < closureFull.size(); iClosure++) {
std::vector<int> &cl = closureFull[iClosure];
for (int iEdge = 0; edges[iEdge] >= 0; iEdge += 2) {
int n0 = cl[edges[iEdge]];
@@ -1202,12 +1203,12 @@ static void generateFaceClosureTetFull(polynomialBasis::clCont &closureFull, std
//Mapping for the p1 nodes
polynomialBasis::clCont closure;
generateFaceClosureTet(closure, 1);
- for (int i = 0; i < closureFull.size(); i++) {
+ for (unsigned int i = 0; i < closureFull.size(); i++) {
std::vector<int> &clFull = closureFull[i];
std::vector<int> &cl = closure[i];
clFull.resize(4, -1);
closureRef[i] = 0;
- for (int j = 0; j < cl.size(); j ++)
+ for (unsigned int j = 0; j < cl.size(); j ++)
clFull[closure[0][j]] = cl[j];
for (int j = 0; j < 4; j ++)
if (clFull[j] == -1)
@@ -1228,7 +1229,7 @@ static void generateFaceClosureTetFull(polynomialBasis::clCont &closureFull, std
if (order >= 3)
generate2dEdgeClosureFull(closureTriangles, closureTrianglesRef, order - 3, 3, false);
addEdgeNodes(closureFull, edges, order);
- for (int iClosure = 0; iClosure < closureFull.size(); iClosure++) {
+ for (unsigned int iClosure = 0; iClosure < closureFull.size(); iClosure++) {
//faces
std::vector<int> &cl = closureFull[iClosure];
if (order >= 3) {
@@ -1241,7 +1242,8 @@ static void generateFaceClosureTetFull(polynomialBasis::clCont &closureFull, std
short int idFace = id/6;
int nOnFace = closureTriangles[iTriClosure].size();
for (int i = 0; i < nOnFace; i++) {
- cl.push_back(4 + 6 * (order - 1) + idFace * nOnFace + closureTriangles[iTriClosure][i]);
+ cl.push_back(4 + 6 * (order - 1) + idFace * nOnFace +
+ closureTriangles[iTriClosure][i]);
}
}
}
@@ -1250,9 +1252,9 @@ static void generateFaceClosureTetFull(polynomialBasis::clCont &closureFull, std
polynomialBasis::clCont insideClosure;
std::vector<int> fakeClosureRef;
generateFaceClosureTetFull(insideClosure, fakeClosureRef, order - 4, false);
- for (int i = 0; i < closureFull.size(); i ++) {
- int shift = closureFull[i].size();
- for (int j = 0; j < insideClosure[i].size(); j++)
+ for (unsigned int i = 0; i < closureFull.size(); i ++) {
+ unsigned int shift = closureFull[i].size();
+ for (unsigned int j = 0; j < insideClosure[i].size(); j++)
closureFull[i].push_back(insideClosure[i][j] + shift);
}
}
@@ -1317,7 +1319,7 @@ static void generateFaceClosurePrismFull(polynomialBasis::clCont &closureFull,
closureRef.resize(40);
generateFaceClosurePrism(closure, 1);
int ref3 = -1, ref4a = -1, ref4b = -1;
- for (int i = 0; i < closure.size(); i++) {
+ for (unsigned int i = 0; i < closure.size(); i++) {
std::vector<int> &clFull = closureFull[i];
std::vector<int> &cl = closure[i];
clFull.resize(6, -1);
@@ -1325,7 +1327,7 @@ static void generateFaceClosurePrismFull(polynomialBasis::clCont &closureFull,
if (ref == -1)
ref = i;
closureRef[i] = ref;
- for (int j = 0; j < cl.size(); j ++)
+ for (unsigned int j = 0; j < cl.size(); j ++)
clFull[closure[ref][j]] = cl[j];
for (int j = 0; j < 6; j ++) {
if (clFull[j] == -1) {
@@ -1335,12 +1337,14 @@ static void generateFaceClosurePrismFull(polynomialBasis::clCont &closureFull,
}
}
}
- static const int edges[] = {0, 1, 0, 2, 0, 3, 1, 2, 1, 4, 2, 5, 3, 4, 3, 5, 4, 5, -1};
+ static const int edges[] = {0, 1, 0, 2, 0, 3, 1, 2, 1, 4, 2, 5,
+ 3, 4, 3, 5, 4, 5, -1};
addEdgeNodes(closureFull, edges, order);
if ( order < 2 )
return;
// face center nodes for p2 prism
- static const int faces[5][4] = {{0, 2, 1, -1}, {3, 4, 5, -1}, {0, 1, 4, 3}, {0, 3, 5, 2}, {1, 2, 5, 4}};
+ static const int faces[5][4] = {{0, 2, 1, -1}, {3, 4, 5, -1}, {0, 1, 4, 3},
+ {0, 3, 5, 2}, {1, 2, 5, 4}};
if ( order == 2 ) {
int nextFaceNode = 15;
@@ -1354,11 +1358,12 @@ static void generateFaceClosurePrismFull(polynomialBasis::clCont &closureFull,
}
nodeSum2Face[nodeSum] = iFace;
}
- for (int i = 0; i < closureFull.size(); i++ ) {
+ for (unsigned int i = 0; i < closureFull.size(); i++ ) {
for (int iFace = 0; iFace < numFaces; iFace++ ) {
int nodeSum = 0;
for (int iNode = 0; iNode < numFaceNodes; iNode++)
- nodeSum += faces[iFace][iNode] < 0 ? faces[iFace][iNode] : closureFull[i][ faces[iFace][iNode] ];
+ nodeSum += faces[iFace][iNode] < 0 ? faces[iFace][iNode] :
+ closureFull[i][ faces[iFace][iNode] ];
std::map<int,int>::iterator it = nodeSum2Face.find(nodeSum);
if (it == nodeSum2Face.end() )
Msg::Error("Prism face not found");
diff --git a/Post/PViewVertexArrays.cpp b/Post/PViewVertexArrays.cpp
index 0ac996ad4b0e65d6c8033c3cbaaa08e89e689d0a..5b31ecbe6aefcaff5b5828aa3a5bb44446efd51b 100644
--- a/Post/PViewVertexArrays.cpp
+++ b/Post/PViewVertexArrays.cpp
@@ -1037,7 +1037,7 @@ static void addTensorElement(PView *p, int iEnt, int iEle, int numNodes, int typ
else if (opt->tensorType == PViewOptions::Ellipse ||
opt->tensorType == PViewOptions::Ellipsoid) {
if(opt->glyphLocation == PViewOptions::Vertex){
- double vval[3][4]= {0,0,0, 0,0,0, 0,0,0, 0,0,0};
+ double vval[3][4]= {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}};
for(int i = 0; i < numNodes; i++){
for (int j = 0; j < 3; j++) {
tensor(j,0) = val [i][0+j*3];
@@ -1060,7 +1060,7 @@ static void addTensorElement(PView *p, int iEnt, int iEle, int numNodes, int typ
}
}
else if(opt->glyphLocation == PViewOptions::COG){
- double vval[3][4]= {0,0,0, 0,0,0, 0,0,0, 0,0,0};
+ double vval[3][4]= {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}};
for(int i = 0; i < numNodes; i++) {
for (int j = 0; j < 3; j++) {
tensor(j,0) = val [i][0+j*3];