diff --git a/Mesh/meshGFaceRecombine.cpp b/Mesh/meshGFaceRecombine.cpp
index 60d8315484c15132362ca8f12d67db63ecd943db..e105c89952951e8226a3345904863df0a59c0c15 100644
--- a/Mesh/meshGFaceRecombine.cpp
+++ b/Mesh/meshGFaceRecombine.cpp
@@ -56,104 +56,107 @@ Rec2DData *Rec2DData::_cur = NULL;
double **Rec2DVertex::_qualVSnum = NULL;
double **Rec2DVertex::_gains = NULL;
+namespace {
+
#ifdef REC2D_DRAW
-static void __draw(double dt = .0)
-{
- double time = Cpu();
- Recombine2D::drawStateOrigin();
- CTX::instance()->mesh.changed = ENT_ALL;
- drawContext::global()->draw();
- while (Cpu()-time < dt)
- FlGui::instance()->check();
-}
+ void __draw(double dt = .0)
+ {
+ double time = Cpu();
+ Recombine2D::drawStateOrigin();
+ CTX::instance()->mesh.changed = ENT_ALL;
+ drawContext::global()->draw();
+ while (Cpu()-time < dt)
+ FlGui::instance()->check();
+ }
-static void __drawWait(double time, double dt)
-{
- Recombine2D::drawStateOrigin();
- CTX::instance()->mesh.changed = ENT_ALL;
- drawContext::global()->draw();
- while (Cpu()-time < dt)
- FlGui::instance()->check();
-}
+ void __drawWait(double time, double dt)
+ {
+ Recombine2D::drawStateOrigin();
+ CTX::instance()->mesh.changed = ENT_ALL;
+ drawContext::global()->draw();
+ while (Cpu()-time < dt)
+ FlGui::instance()->check();
+ }
#endif
-static bool edgesAreInOrder(Rec2DEdge const*const*const edges, const int numEdges)
-{
- Rec2DVertex **v, *v0, *v1;
- v = new Rec2DVertex*[numEdges];
- v0 = edges[0]->getVertex(0);
- v1 = edges[0]->getVertex(1);
- if (edges[1]->getVertex(0) == v0 || edges[1]->getVertex(1) == v0) {
- v[0] = v0;
- if (edges[1]->getVertex(0) == v0)
- v[1] = edges[1]->getVertex(1);
- else
- v[1] = edges[1]->getVertex(0);
- }
- else if (edges[1]->getVertex(0) == v1 || edges[1]->getVertex(1) == v1) {
- v[0] = v1;
- if (edges[1]->getVertex(0) == v1)
- v[1] = edges[1]->getVertex(1);
- else
- v[1] = edges[1]->getVertex(0);
- }
- else {
- Msg::Error("edges not in order (1)");
- for (int i = 0; i < numEdges; ++i) {
- edges[i]->print();
+ bool edgesAreInOrder(Rec2DEdge const*const*const edges, const int numEdges)
+ {
+ Rec2DVertex **v, *v0, *v1;
+ v = new Rec2DVertex*[numEdges];
+ v0 = edges[0]->getVertex(0);
+ v1 = edges[0]->getVertex(1);
+ if (edges[1]->getVertex(0) == v0 || edges[1]->getVertex(1) == v0) {
+ v[0] = v0;
+ if (edges[1]->getVertex(0) == v0)
+ v[1] = edges[1]->getVertex(1);
+ else
+ v[1] = edges[1]->getVertex(0);
+ }
+ else if (edges[1]->getVertex(0) == v1 || edges[1]->getVertex(1) == v1) {
+ v[0] = v1;
+ if (edges[1]->getVertex(0) == v1)
+ v[1] = edges[1]->getVertex(1);
+ else
+ v[1] = edges[1]->getVertex(0);
}
- return false;
- }
- for (int i = 2; i < numEdges; ++i) {
- if (edges[i]->getVertex(0) == v[i-1])
- v[i] = edges[i]->getVertex(1);
- else if (edges[i]->getVertex(1) == v[i-1])
- v[i] = edges[i]->getVertex(0);
else {
- Msg::Error("edges not in order (2)");
+ Msg::Error("edges not in order (1)");
for (int i = 0; i < numEdges; ++i) {
edges[i]->print();
}
return false;
}
- }
- if ((v[0] == v1 && v[numEdges-1] != v0) ||
- (v[0] == v0 && v[numEdges-1] != v1) ) {
- Msg::Error("edges not in order (3)");
- for (int i = 0; i < numEdges; ++i) {
- edges[i]->print();
+ for (int i = 2; i < numEdges; ++i) {
+ if (edges[i]->getVertex(0) == v[i-1])
+ v[i] = edges[i]->getVertex(1);
+ else if (edges[i]->getVertex(1) == v[i-1])
+ v[i] = edges[i]->getVertex(0);
+ else {
+ Msg::Error("edges not in order (2)");
+ for (int i = 0; i < numEdges; ++i) {
+ edges[i]->print();
+ }
+ return false;
+ }
}
- return false;
+ if ((v[0] == v1 && v[numEdges-1] != v0) ||
+ (v[0] == v0 && v[numEdges-1] != v1) ) {
+ Msg::Error("edges not in order (3)");
+ for (int i = 0; i < numEdges; ++i) {
+ edges[i]->print();
+ }
+ return false;
+ }
+ delete v;
+ return true;
}
- delete v;
- return true;
-}
-
-static void __crash()
-{
- Msg::Info(" ");
- Recombine2D::drawStateOrigin();
- int a[2];
- int e = 0;
- for (int i = 0; i < 10000000; ++i) e+=a[i];
- Msg::Info("%d",e);
-}
-//static void __wait(double dt = REC2D_WAIT_TM_3)
-//{
-//#ifdef REC2D_DRAW
-// Msg::Info(" ");
-// double time = Cpu();
-// while (Cpu()-time < dt)
-// FlGui::instance()->check();
-//#endif
-//}
-//
-static int otherParity(const int a)
-{
- if (a % 2)
- return a - 1;
- return a + 1;
+ void __crash()
+ {
+ Msg::Info(" ");
+ Recombine2D::drawStateOrigin();
+ int a[2];
+ int e = 0;
+ for (int i = 0; i < 10000000; ++i) e+=a[i];
+ Msg::Info("%d",e);
+ }
+
+ //void __wait(double dt = REC2D_WAIT_TM_3)
+ //{
+ //#ifdef REC2D_DRAW
+ // Msg::Info(" ");
+ // double time = Cpu();
+ // while (Cpu()-time < dt)
+ // FlGui::instance()->check();
+ //#endif
+ //}
+ //
+ int otherParity(const int a)
+ {
+ if (a % 2)
+ return a - 1;
+ return a + 1;
+ }
}
namespace std // overload of std::swap(..) for Rec2DData::Action class
diff --git a/Numeric/ElementType.h b/Numeric/ElementType.h
index 3338abc6d8e6a0e38e326f05a9fb28107a7b9c83..f7a9073420fc103d5b4b7924473cc7895c0d3864 100644
--- a/Numeric/ElementType.h
+++ b/Numeric/ElementType.h
@@ -6,21 +6,21 @@
#ifndef _ELEMENTTYPE_H_
#define _ELEMENTTYPE_H_
-class ElementType
+namespace ElementType
{
- public :
// Give parent type, order & dimension
// corresponding to any element type.
- static int ParentTypeFromTag(int tag);
- static int OrderFromTag(int tag);
- static int DimensionFromTag(int tag);
+ int ParentTypeFromTag(int tag);
+ int OrderFromTag(int tag);
+ int DimensionFromTag(int tag);
- // Gives > 0 if element type is in Serendipity Family.
- // Gives < 2 if element type is in 'Normal' Family.
+ // Gives > 0 if element tag is in Serendipity Family.
+ // Gives < 2 if element tag is in 'Normal' Family.
// 1 is for element that is either Serendipity or not !
- static int SerendipityFromTag(int tag);
+ int SerendipityFromTag(int tag);
- static int getTag(int parentTag, int order, bool serendip = false);
+ // Give element tag from type, order & serendip
+ int getTag(int parentTag, int order, bool serendip = false);
};
#endif
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index 8975c7006dfcfda5547138548ba55d7db588cebb..9f956eb87c6488dcb2685e0137fd032726c5fa4d 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -13,6 +13,15 @@
#include "BasisFactory.h"
#include "Numeric.h"
+namespace {
+ inline double calcDet3D(double dxdX, double dydX, double dzdX,
+ double dxdY, double dydY, double dzdY,
+ double dxdZ, double dydZ, double dzdZ)
+ {
+ return dxdX*dydY*dzdZ + dxdY*dydZ*dzdX + dydX*dzdY*dxdZ
+ - dxdZ*dydY*dzdX - dxdY*dydX*dzdZ - dydZ*dzdY*dxdX;
+ }
+}
JacobianBasis::JacobianBasis(int tag)
{
@@ -167,14 +176,6 @@ double JacobianBasis::getPrimNormal2D(const fullMatrix<double> &nodesXYZ, fullMa
}
-static inline double calcDet3D(double dxdX, double dydX, double dzdX,
- double dxdY, double dydY, double dzdY,
- double dxdZ, double dydZ, double dzdZ)
-{
- return dxdX*dydY*dzdZ + dxdY*dydZ*dzdX + dydX*dzdY*dxdZ
- - dxdZ*dydY*dzdX - dxdY*dydX*dzdZ - dydZ*dzdY*dxdX;
-}
-
// Returns absolute value of Jacobian of straight volume element at barycenter
double JacobianBasis::getPrimJac3D(const fullMatrix<double> &nodesXYZ) const
{
diff --git a/Numeric/bezierBasis.cpp b/Numeric/bezierBasis.cpp
index cf2d14c8e413d63e9b30f4c83ba82de74c950fc6..b6e60abf36d15465c88611b5d2ff5a989f42f1b7 100644
--- a/Numeric/bezierBasis.cpp
+++ b/Numeric/bezierBasis.cpp
@@ -13,223 +13,184 @@
#include "BasisFactory.h"
#include "Numeric.h"
+namespace {
// Sub Control Points
-static std::vector< fullMatrix<double> > generateSubPointsLine(int order)
-{
- std::vector< fullMatrix<double> > subPoints(2);
-
- subPoints[0] = gmshGenerateMonomialsLine(order);
- subPoints[0].scale(.5/order);
-
- 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] = gmshGenerateMonomialsTriangle(order);
- subPoints[0].scale(.5/order);
+ std::vector< fullMatrix<double> > generateSubPointsLine(int order)
+ {
+ std::vector< fullMatrix<double> > subPoints(2);
- subPoints[1].copy(subPoints[0]);
- prox.setAsProxy(subPoints[1], 0, 1);
- prox.add(.5);
+ subPoints[0] = gmshGenerateMonomialsLine(order);
+ subPoints[0].scale(.5/order);
- subPoints[2].copy(subPoints[0]);
- prox.setAsProxy(subPoints[2], 1, 1);
- prox.add(.5);
+ subPoints[1].copy(subPoints[0]);
+ subPoints[1].add(.5);
- subPoints[3].copy(subPoints[0]);
- subPoints[3].scale(-1.);
- subPoints[3].add(.5);
+ return subPoints;
+ }
- return subPoints;
-}
+ std::vector< fullMatrix<double> > generateSubPointsTriangle(int order)
+ {
+ std::vector< fullMatrix<double> > subPoints(4);
+ fullMatrix<double> prox;
-static std::vector< fullMatrix<double> > generateSubPointsQuad(int order)
-{
- std::vector< fullMatrix<double> > subPoints(4);
- fullMatrix<double> prox;
+ subPoints[0] = gmshGenerateMonomialsTriangle(order);
+ subPoints[0].scale(.5/order);
- subPoints[0] = gmshGenerateMonomialsQuadrangle(order);
- subPoints[0].scale(.5/order);
+ subPoints[1].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[1], 0, 1);
+ prox.add(.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[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);
- subPoints[3].copy(subPoints[1]);
- prox.setAsProxy(subPoints[3], 1, 1);
- prox.add(.5);
+ return subPoints;
+ }
- return subPoints;
-}
+ std::vector< fullMatrix<double> > generateSubPointsQuad(int order)
+ {
+ std::vector< fullMatrix<double> > subPoints(4);
+ fullMatrix<double> prox;
-static std::vector< fullMatrix<double> > generateSubPointsTetrahedron(int order)
-{
- std::vector< fullMatrix<double> > subPoints(8);
- fullMatrix<double> prox1;
- fullMatrix<double> prox2;
-
- subPoints[0] = gmshGenerateMonomialsTetrahedron(order);
- subPoints[0].scale(.5/order);
-
- 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;
-}
+ subPoints[0] = gmshGenerateMonomialsQuadrangle(order);
+ subPoints[0].scale(.5/order);
-static std::vector< fullMatrix<double> > generateSubPointsPrism(int order)
-{
- std::vector< fullMatrix<double> > subPoints(8);
- fullMatrix<double> prox;
+ subPoints[1].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[1], 0, 1);
+ prox.add(.5);
- subPoints[0] = gmshGenerateMonomialsPrism(order);
- subPoints[0].scale(.5/order);
+ subPoints[2].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[2], 1, 1);
+ prox.add(.5);
- subPoints[1].copy(subPoints[0]);
- prox.setAsProxy(subPoints[1], 0, 1);
- prox.add(.5);
+ subPoints[3].copy(subPoints[1]);
+ prox.setAsProxy(subPoints[3], 1, 1);
+ prox.add(.5);
- subPoints[2].copy(subPoints[0]);
- prox.setAsProxy(subPoints[2], 1, 1);
- prox.add(.5);
+ return subPoints;
+ }
- subPoints[3].copy(subPoints[0]);
- prox.setAsProxy(subPoints[3], 0, 2);
- prox.scale(-1.);
- prox.add(.5);
+ std::vector< fullMatrix<double> > generateSubPointsTetrahedron(int order)
+ {
+ std::vector< fullMatrix<double> > subPoints(8);
+ fullMatrix<double> prox1;
+ fullMatrix<double> prox2;
- subPoints[4].copy(subPoints[0]);
- prox.setAsProxy(subPoints[4], 2, 1);
- prox.add(.5);
+ subPoints[0] = gmshGenerateMonomialsTetrahedron(order);
+ subPoints[0].scale(.5/order);
- subPoints[5].copy(subPoints[1]);
- prox.setAsProxy(subPoints[5], 2, 1);
- prox.add(.5);
+ subPoints[1].copy(subPoints[0]);
+ prox1.setAsProxy(subPoints[1], 0, 1);
+ prox1.add(.5);
- subPoints[6].copy(subPoints[2]);
- prox.setAsProxy(subPoints[6], 2, 1);
- prox.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);
- subPoints[7].copy(subPoints[3]);
- prox.setAsProxy(subPoints[7], 2, 1);
- prox.add(.5);
- return subPoints;
-}
+ return subPoints;
+ }
-static std::vector< fullMatrix<double> > generateSubPointsHex(int order)
-{
- std::vector< fullMatrix<double> > subPoints(8);
- fullMatrix<double> prox;
+ std::vector< fullMatrix<double> > generateSubPointsPrism(int order)
+ {
+ std::vector< fullMatrix<double> > subPoints(8);
+ fullMatrix<double> prox;
- subPoints[0] = gmshGenerateMonomialsHexahedron(order);
- subPoints[0].scale(.5/order);
+ subPoints[0] = gmshGenerateMonomialsPrism(order);
+ subPoints[0].scale(.5/order);
- subPoints[1].copy(subPoints[0]);
- prox.setAsProxy(subPoints[1], 0, 1);
- prox.add(.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[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);
+ 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[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[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[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);
+ subPoints[7].copy(subPoints[3]);
+ prox.setAsProxy(subPoints[7], 2, 1);
+ prox.add(.5);
- return subPoints;
-}
+ return subPoints;
+ }
-static std::vector< fullMatrix<double> > generateSubPointsPyr(int order)
-{
- if (order == 0) {
- std::vector< fullMatrix<double> > subPoints(4);
+ std::vector< fullMatrix<double> > generateSubPointsHex(int order)
+ {
+ std::vector< fullMatrix<double> > subPoints(8);
fullMatrix<double> prox;
- subPoints[0] = JacobianBasis::generateJacMonomialsPyramid(order);
- subPoints[0].scale(.5/(order+2));
+ subPoints[0] = gmshGenerateMonomialsHexahedron(order);
+ subPoints[0].scale(.5/order);
subPoints[1].copy(subPoints[0]);
prox.setAsProxy(subPoints[1], 0, 1);
@@ -243,207 +204,248 @@ static std::vector< fullMatrix<double> > generateSubPointsPyr(int order)
prox.setAsProxy(subPoints[3], 1, 1);
prox.add(.5);
- return subPoints;
- }
+ 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);
- std::vector< fullMatrix<double> > subPoints(8);
- fullMatrix<double> ref, prox;
+ subPoints[7].copy(subPoints[3]);
+ prox.setAsProxy(subPoints[7], 2, 1);
+ prox.add(.5);
- subPoints[0] = JacobianBasis::generateJacMonomialsPyramid(order);
- int nPts = subPoints[0].size1();
- for (int i = 0; i < nPts; ++i) {
- const double factor = .5 / (order + 2 - subPoints[0](i, 2));
- subPoints[0](i, 0) = subPoints[0](i, 0) * factor;
- subPoints[0](i, 1) = subPoints[0](i, 1) * factor;
- subPoints[0](i, 2) = subPoints[0](i, 2) * .5 / (order + 2);
+ return subPoints;
}
- subPoints[1].copy(subPoints[0]);
- prox.setAsProxy(subPoints[1], 0, 1);
- prox.add(.5);
+ std::vector< fullMatrix<double> > generateSubPointsPyr(int order)
+ {
+ if (order == 0) {
+ std::vector< fullMatrix<double> > subPoints(4);
+ fullMatrix<double> prox;
- subPoints[2].copy(subPoints[0]);
- prox.setAsProxy(subPoints[2], 1, 1);
- prox.add(.5);
+ subPoints[0] = JacobianBasis::generateJacMonomialsPyramid(order);
+ subPoints[0].scale(.5/(order+2));
- subPoints[3].copy(subPoints[1]);
- prox.setAsProxy(subPoints[3], 1, 1);
- prox.add(.5);
+ subPoints[1].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[1], 0, 1);
+ prox.add(.5);
- subPoints[4].copy(subPoints[0]);
- prox.setAsProxy(subPoints[4], 2, 1);
- prox.add(.5);
+ subPoints[2].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[2], 1, 1);
+ prox.add(.5);
- subPoints[5].copy(subPoints[1]);
- prox.setAsProxy(subPoints[5], 2, 1);
- prox.add(.5);
+ subPoints[3].copy(subPoints[1]);
+ prox.setAsProxy(subPoints[3], 1, 1);
+ prox.add(.5);
- subPoints[6].copy(subPoints[2]);
- prox.setAsProxy(subPoints[6], 2, 1);
- prox.add(.5);
+ return subPoints;
+ }
- subPoints[7].copy(subPoints[3]);
- prox.setAsProxy(subPoints[7], 2, 1);
- prox.add(.5);
+ std::vector< fullMatrix<double> > subPoints(8);
+ fullMatrix<double> ref, prox;
- for (int i = 0; i < 8; ++i) {
- for (int j = 0; j < nPts; ++j) {
- const double factor = (1. - subPoints[i](j, 2));
- subPoints[i](j, 0) = subPoints[i](j, 0) * factor;
- subPoints[i](j, 1) = subPoints[i](j, 1) * factor;
+ subPoints[0] = JacobianBasis::generateJacMonomialsPyramid(order);
+ int nPts = subPoints[0].size1();
+ for (int i = 0; i < nPts; ++i) {
+ const double factor = .5 / (order + 2 - subPoints[0](i, 2));
+ subPoints[0](i, 0) = subPoints[0](i, 0) * factor;
+ subPoints[0](i, 1) = subPoints[0](i, 1) * factor;
+ subPoints[0](i, 2) = subPoints[0](i, 2) * .5 / (order + 2);
}
- }
- return subPoints;
-}
+ subPoints[1].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[1], 0, 1);
+ prox.add(.5);
-// Matrices generation
-static int nChoosek(int n, int k)
-{
- if (n < k || k < 0) {
- Msg::Error("Wrong argument for combination.");
- return 1;
+ 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);
+
+ 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);
+
+ for (int i = 0; i < 8; ++i) {
+ for (int j = 0; j < nPts; ++j) {
+ const double factor = (1. - subPoints[i](j, 2));
+ subPoints[i](j, 0) = subPoints[i](j, 0) * factor;
+ subPoints[i](j, 1) = subPoints[i](j, 1) * factor;
+ }
+ }
+
+ return subPoints;
}
- if (k > n/2) k = n-k;
- if (k == 1)
- return n;
- if (k == 0)
- return 1;
+ // Matrices generation
+ int nChoosek(int n, int k)
+ {
+ if (n < k || k < 0) {
+ Msg::Error("Wrong argument for combination.");
+ return 1;
+ }
- int c = 1;
- for (int i = 1; i <= k; i++, n--) (c *= n) /= i;
- return c;
-}
+ if (k > n/2) k = n-k;
+ if (k == 1)
+ return n;
+ if (k == 0)
+ return 1;
-static fullMatrix<double> generateBez2LagMatrix
- (const fullMatrix<double> &exponent, const fullMatrix<double> &point,
- int order, int dimSimplex)
-{
- if(exponent.size1() != point.size1() || exponent.size2() != point.size2()){
- Msg::Fatal("Wrong sizes for bez2lag matrix generation %d %d -- %d %d",
- exponent.size1(),point.size1(),
- exponent.size2(),point.size2());
- return fullMatrix<double>(1, 1);
+ int c = 1;
+ for (int i = 1; i <= k; i++, n--) (c *= n) /= i;
+ return c;
}
- int ndofs = exponent.size1();
- int dim = exponent.size2();
-
- fullMatrix<double> bez2Lag(ndofs, ndofs);
- for (int i = 0; i < ndofs; i++) {
- for (int j = 0; j < ndofs; j++) {
- double dd = 1.;
-
- {
- double pointCompl = 1.;
- int exponentCompl = order;
- for (int k = 0; k < dimSimplex; k++) {
- dd *= nChoosek(exponentCompl, (int) exponent(i, k))
- * pow(point(j, k), exponent(i, k));
- pointCompl -= point(j, k);
- exponentCompl -= (int) exponent(i, k);
+ fullMatrix<double> generateBez2LagMatrix
+ (const fullMatrix<double> &exponent, const fullMatrix<double> &point,
+ int order, int dimSimplex)
+ {
+ if(exponent.size1() != point.size1() || exponent.size2() != point.size2()){
+ Msg::Fatal("Wrong sizes for bez2lag matrix generation %d %d -- %d %d",
+ exponent.size1(),point.size1(),
+ exponent.size2(),point.size2());
+ return fullMatrix<double>(1, 1);
+ }
+
+ int ndofs = exponent.size1();
+ int dim = exponent.size2();
+
+ fullMatrix<double> bez2Lag(ndofs, ndofs);
+ for (int i = 0; i < ndofs; i++) {
+ for (int j = 0; j < ndofs; j++) {
+ double dd = 1.;
+
+ {
+ double pointCompl = 1.;
+ int exponentCompl = order;
+ for (int k = 0; k < dimSimplex; k++) {
+ dd *= nChoosek(exponentCompl, (int) exponent(i, k))
+ * pow(point(j, k), exponent(i, k));
+ pointCompl -= point(j, k);
+ exponentCompl -= (int) exponent(i, k);
+ }
+ dd *= pow(pointCompl, exponentCompl);
}
- dd *= pow(pointCompl, exponentCompl);
- }
- for (int k = dimSimplex; k < dim; k++)
- dd *= nChoosek(order, (int) exponent(i, k))
- * pow(point(j, k), exponent(i, k))
- * pow(1. - point(j, k), order - exponent(i, k));
+ for (int k = dimSimplex; k < dim; k++)
+ dd *= nChoosek(order, (int) exponent(i, k))
+ * pow(point(j, k), exponent(i, k))
+ * pow(1. - point(j, k), order - exponent(i, k));
- bez2Lag(j, i) = dd;
+ bez2Lag(j, i) = dd;
+ }
}
+ return bez2Lag;
}
- return bez2Lag;
-}
-static fullMatrix<double> generateBez2LagMatrixPyramid
- (const fullMatrix<double> &exponent, const fullMatrix<double> &point,
- int order)
-{
- if(exponent.size1() != point.size1() || exponent.size2() != point.size2() ||
- exponent.size2() != 3){
- Msg::Fatal("Wrong sizes for bez2lag matrix generation %d %d -- %d %d",
- exponent.size1(),point.size1(),
- exponent.size2(),point.size2());
- return fullMatrix<double>(1, 1);
- }
+ fullMatrix<double> generateBez2LagMatrixPyramid
+ (const fullMatrix<double> &exponent, const fullMatrix<double> &point,
+ int order)
+ {
+ if(exponent.size1() != point.size1() || exponent.size2() != point.size2() ||
+ exponent.size2() != 3){
+ Msg::Fatal("Wrong sizes for bez2lag matrix generation %d %d -- %d %d",
+ exponent.size1(),point.size1(),
+ exponent.size2(),point.size2());
+ return fullMatrix<double>(1, 1);
+ }
- int ndofs = exponent.size1();
-
- fullMatrix<double> bez2Lag(ndofs, ndofs);
- for (int i = 0; i < ndofs; i++) {
- for (int j = 0; j < ndofs; j++) {
- bez2Lag(i, j) =
- nChoosek(order, exponent(j, 2))
- * pow(point(i, 2), exponent(j, 2)) //FIXME use pow_int
- * pow(1. - point(i, 2), order - exponent(j, 2));
-
- double p = order + 2 - exponent(j, 2);
- double denom = 1. - point(i, 2);
- bez2Lag(i, j) *=
- nChoosek(p, exponent(j, 0))
- * nChoosek(p, exponent(j, 1))
- * pow(point(i, 0) / denom, exponent(j, 0))
- * pow(point(i, 1) / denom, exponent(j, 1))
- * pow(1. - point(i, 0) / denom, p - exponent(j, 0))
- * pow(1. - point(i, 1) / denom, p - exponent(j, 1));
+ int ndofs = exponent.size1();
+
+ fullMatrix<double> bez2Lag(ndofs, ndofs);
+ for (int i = 0; i < ndofs; i++) {
+ for (int j = 0; j < ndofs; j++) {
+ bez2Lag(i, j) =
+ nChoosek(order, exponent(j, 2))
+ * pow(point(i, 2), exponent(j, 2)) //FIXME use pow_int
+ * pow(1. - point(i, 2), order - exponent(j, 2));
+
+ double p = order + 2 - exponent(j, 2);
+ double denom = 1. - point(i, 2);
+ bez2Lag(i, j) *=
+ nChoosek(p, exponent(j, 0))
+ * nChoosek(p, exponent(j, 1))
+ * pow(point(i, 0) / denom, exponent(j, 0))
+ * pow(point(i, 1) / denom, exponent(j, 1))
+ * pow(1. - point(i, 0) / denom, p - exponent(j, 0))
+ * pow(1. - point(i, 1) / denom, p - exponent(j, 1));
+ }
}
+ return bez2Lag;
}
- return bez2Lag;
-}
-static fullMatrix<double> generateSubDivisor
- (const fullMatrix<double> &exponents, const std::vector< fullMatrix<double> > &subPoints,
- const fullMatrix<double> &lag2Bez, int order, int dimSimplex)
-{
- if (exponents.size1() != lag2Bez.size1() || exponents.size1() != lag2Bez.size2()){
- Msg::Fatal("Wrong sizes for Bezier Divisor %d %d -- %d %d",
- exponents.size1(), lag2Bez.size1(),
- exponents.size1(), lag2Bez.size2());
- return fullMatrix<double>(1, 1);
- }
+ fullMatrix<double> generateSubDivisor
+ (const fullMatrix<double> &exponents, const std::vector< fullMatrix<double> > &subPoints,
+ const fullMatrix<double> &lag2Bez, int order, int dimSimplex)
+ {
+ if (exponents.size1() != lag2Bez.size1() || exponents.size1() != lag2Bez.size2()){
+ Msg::Fatal("Wrong sizes for Bezier Divisor %d %d -- %d %d",
+ exponents.size1(), lag2Bez.size1(),
+ exponents.size1(), lag2Bez.size2());
+ return fullMatrix<double>(1, 1);
+ }
- int nbPts = lag2Bez.size1();
- int nbSubPts = nbPts * subPoints.size();
+ int nbPts = lag2Bez.size1();
+ int nbSubPts = nbPts * subPoints.size();
- fullMatrix<double> intermediate2(nbPts, nbPts);
- fullMatrix<double> subDivisor(nbSubPts, nbPts);
+ fullMatrix<double> intermediate2(nbPts, nbPts);
+ fullMatrix<double> subDivisor(nbSubPts, nbPts);
- for (unsigned int i = 0; i < subPoints.size(); i++) {
- fullMatrix<double> intermediate1 =
- generateBez2LagMatrix(exponents, subPoints[i], order, dimSimplex);
- lag2Bez.mult(intermediate1, intermediate2);
- subDivisor.copy(intermediate2, 0, nbPts, 0, nbPts, i*nbPts, 0);
+ for (unsigned int i = 0; i < subPoints.size(); i++) {
+ fullMatrix<double> intermediate1 =
+ generateBez2LagMatrix(exponents, subPoints[i], order, dimSimplex);
+ lag2Bez.mult(intermediate1, intermediate2);
+ subDivisor.copy(intermediate2, 0, nbPts, 0, nbPts, i*nbPts, 0);
+ }
+ return subDivisor;
}
- return subDivisor;
-}
-static fullMatrix<double> generateSubDivisorPyramid
- (const fullMatrix<double> &exponents, const std::vector< fullMatrix<double> > &subPoints,
- const fullMatrix<double> &lag2Bez, int order)
-{
- if (exponents.size1() != lag2Bez.size1() || exponents.size1() != lag2Bez.size2()){
- Msg::Fatal("Wrong sizes for Bezier Divisor %d %d -- %d %d",
- exponents.size1(), lag2Bez.size1(),
- exponents.size1(), lag2Bez.size2());
- return fullMatrix<double>(1, 1);
- }
+ fullMatrix<double> generateSubDivisorPyramid
+ (const fullMatrix<double> &exponents, const std::vector< fullMatrix<double> > &subPoints,
+ const fullMatrix<double> &lag2Bez, int order)
+ {
+ if (exponents.size1() != lag2Bez.size1() || exponents.size1() != lag2Bez.size2()){
+ Msg::Fatal("Wrong sizes for Bezier Divisor %d %d -- %d %d",
+ exponents.size1(), lag2Bez.size1(),
+ exponents.size1(), lag2Bez.size2());
+ return fullMatrix<double>(1, 1);
+ }
- int nbPts = lag2Bez.size1();
- int nbSubPts = nbPts * subPoints.size();
+ int nbPts = lag2Bez.size1();
+ int nbSubPts = nbPts * subPoints.size();
- fullMatrix<double> intermediate2(nbPts, nbPts);
- fullMatrix<double> subDivisor(nbSubPts, nbPts);
+ fullMatrix<double> intermediate2(nbPts, nbPts);
+ fullMatrix<double> subDivisor(nbSubPts, nbPts);
- for (unsigned int i = 0; i < subPoints.size(); i++) {
- fullMatrix<double> intermediate1 =
- generateBez2LagMatrixPyramid(exponents, subPoints[i], order);
- lag2Bez.mult(intermediate1, intermediate2);
- subDivisor.copy(intermediate2, 0, nbPts, 0, nbPts, i*nbPts, 0);
+ for (unsigned int i = 0; i < subPoints.size(); i++) {
+ fullMatrix<double> intermediate1 =
+ generateBez2LagMatrixPyramid(exponents, subPoints[i], order);
+ lag2Bez.mult(intermediate1, intermediate2);
+ subDivisor.copy(intermediate2, 0, nbPts, 0, nbPts, i*nbPts, 0);
+ }
+ return subDivisor;
}
- return subDivisor;
}
void bezierBasis::_construct(int parentType, int p)
diff --git a/Numeric/nodalBasis.cpp b/Numeric/nodalBasis.cpp
index c6c97ea1e936c0609a11636483572e21d819f398..6a8cb4e14100facefb212b77647664e6b1a81a2a 100644
--- a/Numeric/nodalBasis.cpp
+++ b/Numeric/nodalBasis.cpp
@@ -10,591 +10,561 @@
#include "pointsGenerators.h"
-static void generate1dVertexClosure(nodalBasis::clCont &closure, int order)
-{
- closure.clear();
- closure.resize(2);
- closure[0].push_back(0);
- closure[1].push_back(order == 0 ? 0 : 1);
- closure[0].type = MSH_PNT;
- closure[1].type = MSH_PNT;
-}
-
-
-static void generate1dVertexClosureFull(nodalBasis::clCont &closure,
- std::vector<int> &closureRef, int order)
-{
- closure.clear();
- closure.resize(2);
- closure[0].push_back(0);
- if (order != 0) {
- closure[0].push_back(1);
- closure[1].push_back(1);
+namespace ClosureGen {
+ inline double pow2(double x)
+ {
+ return x * x;
}
- closure[1].push_back(0);
- for (int i = 0; i < order - 1; i++) {
- closure[0].push_back(2 + i);
- closure[1].push_back(2 + order - 2 - i);
+
+ void rotateHex(int iFace, int iRot, int iSign, double uI, double vI,
+ double &uO, double &vO, double &wO)
+ {
+ if (iSign<0){
+ double tmp=uI;
+ uI=vI;
+ vI=tmp;
+ }
+ for (int i=0; i < iRot; i++){
+ double tmp=uI;
+ uI=-vI;
+ vI=tmp;
+ }
+ switch (iFace) {
+ case 0: uO = vI; vO = uI; wO=-1; break;
+ case 1: uO = uI; vO = -1; wO=vI; break;
+ case 2: uO =-1; vO = vI; wO=uI; break;
+ case 3: uO = 1; vO = uI; wO=vI; break;
+ case 4: uO =-uI; vO = 1; wO=vI; break;
+ case 5: uO =uI; vO = vI; wO=1; break;
+ }
}
- closureRef.resize(2);
- closureRef[0] = 0;
- closureRef[1] = 0;
-}
+ void rotateHexFull(int iFace, int iRot, int iSign, double uI, double vI, double wI, double &uO, double &vO, double &wO)
+ {
+ switch (iFace) {
+ case 0: uO = uI; vO = vI; wO = wI; break;
+ case 1: uO = wI; vO = uI; wO = vI; break;
+ case 2: uO = vI; vO = wI; wO = uI; break;
+ case 3: uO = wI; vO = vI; wO =-uI; break;
+ case 4: uO = wI; vO =-uI; wO =-vI; break;
+ case 5: uO = vI; vO = uI; wO =-wI; break;
+ }
+ for (int i = 0; i < iRot; i++){
+ double tmp = uO;
+ uO = -vO;
+ vO = tmp;
+ }
+ if (iSign<0){
+ double tmp = uO;
+ uO = vO;
+ vO = tmp;
+ }
+ }
+ void generate1dVertexClosure(nodalBasis::clCont &closure, int order)
+ {
+ closure.clear();
+ closure.resize(2);
+ closure[0].push_back(0);
+ closure[1].push_back(order == 0 ? 0 : 1);
+ closure[0].type = MSH_PNT;
+ closure[1].type = MSH_PNT;
+ }
+
+ void generate1dVertexClosureFull(nodalBasis::clCont &closure,
+ std::vector<int> &closureRef, int order)
+ {
+ closure.clear();
+ closure.resize(2);
+ closure[0].push_back(0);
+ if (order != 0) {
+ closure[0].push_back(1);
+ closure[1].push_back(1);
+ }
+ closure[1].push_back(0);
+ for (int i = 0; i < order - 1; i++) {
+ closure[0].push_back(2 + i);
+ closure[1].push_back(2 + order - 2 - i);
+ }
+ closureRef.resize(2);
+ closureRef[0] = 0;
+ closureRef[1] = 0;
+ }
-static void getFaceClosureTet(int iFace, int iSign, int iRotate,
- nodalBasis::closure &closure, int order)
-{
- closure.clear();
- closure.resize((order + 1) * (order + 2) / 2);
- closure.type = ElementType::getTag(TYPE_TRI, order, false);
+ void getFaceClosureTet(int iFace, int iSign, int iRotate,
+ nodalBasis::closure &closure, int order)
+ {
+ closure.clear();
+ closure.resize((order + 1) * (order + 2) / 2);
+ closure.type = ElementType::getTag(TYPE_TRI, order, false);
- switch (order){
- case 0:
- closure[0] = 0;
- break;
- default:
- int face[4][3] = {{-3, -2, -1}, {1, -6, 4}, {-4, 5, 3}, {6, 2, -5}};
- int order1node[4][3] = {{0, 2, 1}, {0, 1, 3}, {0, 3, 2}, {3, 1, 2}};
- for (int i = 0; i < 3; ++i){
- int k = (3 + (iSign * i) + iRotate) % 3; //- iSign * iRotate
- closure[i] = order1node[iFace][k];
- }
- for (int i = 0;i < 3; ++i){
- int edgenumber = iSign *
- face[iFace][(6 + i * iSign + (-1 + iSign) / 2 + iRotate) % 3]; //- iSign * iRotate
- for (int k = 0; k < (order - 1); k++){
- if (edgenumber > 0)
- closure[3 + i * (order - 1) + k] =
- 4 + (edgenumber - 1) * (order - 1) + k;
- else
- closure[3 + i * (order - 1) + k] =
- 4 + (-edgenumber) * (order - 1) - 1 - k;
+ switch (order){
+ case 0:
+ closure[0] = 0;
+ break;
+ default:
+ int face[4][3] = {{-3, -2, -1}, {1, -6, 4}, {-4, 5, 3}, {6, 2, -5}};
+ int order1node[4][3] = {{0, 2, 1}, {0, 1, 3}, {0, 3, 2}, {3, 1, 2}};
+ for (int i = 0; i < 3; ++i){
+ int k = (3 + (iSign * i) + iRotate) % 3; //- iSign * iRotate
+ closure[i] = order1node[iFace][k];
}
- }
- int fi = 3 + 3 * (order - 1);
- int ti = 4 + 6 * (order - 1);
- int ndofff = (order - 3 + 2) * (order - 3 + 1) / 2;
- ti = ti + iFace * ndofff;
- for (int k = 0; k < order / 3; k++){
- int orderint = order - 3 - k * 3;
- if (orderint > 0){
- for (int ci = 0; ci < 3 ; ci++){
- int shift = (3 + iSign * ci + iRotate) % 3; //- iSign * iRotate
- closure[fi + ci] = ti + shift;
+ for (int i = 0;i < 3; ++i){
+ int edgenumber = iSign *
+ face[iFace][(6 + i * iSign + (-1 + iSign) / 2 + iRotate) % 3]; //- iSign * iRotate
+ for (int k = 0; k < (order - 1); k++){
+ if (edgenumber > 0)
+ closure[3 + i * (order - 1) + k] =
+ 4 + (edgenumber - 1) * (order - 1) + k;
+ else
+ closure[3 + i * (order - 1) + k] =
+ 4 + (-edgenumber) * (order - 1) - 1 - k;
}
- fi = fi + 3; ti = ti + 3;
- for (int l = 0; l < orderint - 1; l++){
- for (int ei = 0; ei < 3; ei++){
- int edgenumber = (6 + ei * iSign + (-1 + iSign) / 2 + iRotate) % 3;
- //- iSign * iRotate
- if (iSign > 0)
- closure[fi + ei * (orderint - 1) + l] =
- ti + edgenumber * (orderint - 1) + l;
- else
- closure[fi + ei * (orderint - 1) + l] =
- ti + (1 + edgenumber) * (orderint - 1) - 1 - l;
+ }
+ int fi = 3 + 3 * (order - 1);
+ int ti = 4 + 6 * (order - 1);
+ int ndofff = (order - 3 + 2) * (order - 3 + 1) / 2;
+ ti = ti + iFace * ndofff;
+ for (int k = 0; k < order / 3; k++){
+ int orderint = order - 3 - k * 3;
+ if (orderint > 0){
+ for (int ci = 0; ci < 3 ; ci++){
+ int shift = (3 + iSign * ci + iRotate) % 3; //- iSign * iRotate
+ closure[fi + ci] = ti + shift;
}
+ fi = fi + 3; ti = ti + 3;
+ for (int l = 0; l < orderint - 1; l++){
+ for (int ei = 0; ei < 3; ei++){
+ int edgenumber = (6 + ei * iSign + (-1 + iSign) / 2 + iRotate) % 3;
+ //- iSign * iRotate
+ if (iSign > 0)
+ closure[fi + ei * (orderint - 1) + l] =
+ ti + edgenumber * (orderint - 1) + l;
+ else
+ closure[fi + ei * (orderint - 1) + l] =
+ ti + (1 + edgenumber) * (orderint - 1) - 1 - l;
+ }
+ }
+ fi = fi + 3 * (orderint - 1); ti = ti + 3 * (orderint - 1);
+ }
+ else {
+ closure[fi] = ti;
+ ti++;
+ fi++;
}
- fi = fi + 3 * (orderint - 1); ti = ti + 3 * (orderint - 1);
- }
- else {
- closure[fi] = ti;
- ti++;
- fi++;
}
+ break;
}
- break;
}
-}
-
-
-static void generate2dEdgeClosureFull(nodalBasis::clCont &closure,
- std::vector<int> &closureRef,
- int order, int nNod, bool serendip)
-{
- closure.clear();
- closure.resize(2*nNod);
- closureRef.resize(2*nNod);
- int shift = 0;
- for (int corder = order; corder>=0; corder -= (nNod == 3 ? 3 : 2)) {
- if (corder == 0) {
+ void generate2dEdgeClosureFull(nodalBasis::clCont &closure,
+ std::vector<int> &closureRef,
+ int order, int nNod, bool serendip)
+ {
+ closure.clear();
+ closure.resize(2*nNod);
+ closureRef.resize(2*nNod);
+ int shift = 0;
+ for (int corder = order; corder>=0; corder -= (nNod == 3 ? 3 : 2)) {
+ if (corder == 0) {
+ for (int r = 0; r < nNod ; r++){
+ closure[r].push_back(shift);
+ closure[r+nNod].push_back(shift);
+ }
+ break;
+ }
for (int r = 0; r < nNod ; r++){
- closure[r].push_back(shift);
- closure[r+nNod].push_back(shift);
+ for (int j = 0; j < nNod; j++) {
+ closure[r].push_back(shift + (r + j) % nNod);
+ closure[r + nNod].push_back(shift + (r - j + 1 + nNod) % nNod);
+ }
}
- break;
- }
- for (int r = 0; r < nNod ; r++){
- for (int j = 0; j < nNod; j++) {
- closure[r].push_back(shift + (r + j) % nNod);
- closure[r + nNod].push_back(shift + (r - j + 1 + nNod) % nNod);
+ shift += nNod;
+ int n = nNod*(corder-1);
+ for (int r = 0; r < nNod ; r++){
+ for (int j = 0; j < n; j++){
+ closure[r].push_back(shift + (j + (corder - 1) * r) % n);
+ closure[r + nNod].push_back(shift + (n - j - 1 + (corder - 1) * (r + 1)) % n);
+ }
}
+ shift += n;
+ if (serendip) break;
}
- shift += nNod;
- int n = nNod*(corder-1);
- for (int r = 0; r < nNod ; r++){
- for (int j = 0; j < n; j++){
- closure[r].push_back(shift + (j + (corder - 1) * r) % n);
- closure[r + nNod].push_back(shift + (n - j - 1 + (corder - 1) * (r + 1)) % n);
- }
+ for (int r = 0; r < nNod*2 ; r++) {
+ closure[r].type = ElementType::getTag(TYPE_LIN, order);
+ closureRef[r] = 0;
}
- shift += n;
- if (serendip) break;
}
- for (int r = 0; r < nNod*2 ; r++) {
- closure[r].type = ElementType::getTag(TYPE_LIN, order);
- closureRef[r] = 0;
- }
-}
-
-
-static void addEdgeNodes(nodalBasis::clCont &closureFull, const int *edges, int order)
-{
- if (order < 2)
- return;
- int numNodes = 0;
- for (int i = 0; edges[i] >= 0; ++i) {
- numNodes = std::max(numNodes, edges[i] + 1);
- }
- std::vector<std::vector<int> > nodes2edges(numNodes, std::vector<int>(numNodes, -1));
- for (int i = 0; edges[i] >= 0; i += 2) {
- nodes2edges[edges[i]][edges[i + 1]] = i;
- nodes2edges[edges[i + 1]][edges[i]] = i + 1;
- }
- for (unsigned int iClosure = 0; iClosure < closureFull.size(); iClosure++) {
- std::vector<int> &cl = closureFull[iClosure];
- for (int iEdge = 0; edges[iEdge] >= 0; iEdge += 2) {
- if (cl.empty())
- continue;
- int n0 = cl[edges[iEdge]];
- int n1 = cl[edges[iEdge + 1]];
- int oEdge = nodes2edges[n0][n1];
- if (oEdge == -1)
- Msg::Error("invalid p1 closure or invalid edges list");
- for (int i = 0 ; i < order - 1; i++)
- cl.push_back(numNodes + oEdge/2 * (order - 1) + ((oEdge % 2) ? order - 2 - i : i));
+ void addEdgeNodes(nodalBasis::clCont &closureFull, const int *edges, int order)
+ {
+ if (order < 2)
+ return;
+ int numNodes = 0;
+ for (int i = 0; edges[i] >= 0; ++i) {
+ numNodes = std::max(numNodes, edges[i] + 1);
+ }
+ std::vector<std::vector<int> > nodes2edges(numNodes, std::vector<int>(numNodes, -1));
+ for (int i = 0; edges[i] >= 0; i += 2) {
+ nodes2edges[edges[i]][edges[i + 1]] = i;
+ nodes2edges[edges[i + 1]][edges[i]] = i + 1;
+ }
+ for (unsigned int iClosure = 0; iClosure < closureFull.size(); iClosure++) {
+ std::vector<int> &cl = closureFull[iClosure];
+ for (int iEdge = 0; edges[iEdge] >= 0; iEdge += 2) {
+ if (cl.empty())
+ continue;
+ int n0 = cl[edges[iEdge]];
+ int n1 = cl[edges[iEdge + 1]];
+ int oEdge = nodes2edges[n0][n1];
+ if (oEdge == -1)
+ Msg::Error("invalid p1 closure or invalid edges list");
+ for (int i = 0 ; i < order - 1; i++)
+ cl.push_back(numNodes + oEdge/2 * (order - 1) + ((oEdge % 2) ? order - 2 - i : i));
+ }
}
}
-}
-
-
-static void generateFaceClosureTet(nodalBasis::clCont &closure, int order)
-{
- closure.clear();
- for (int iRotate = 0; iRotate < 3; iRotate++){
- for (int iSign = 1; iSign >= -1; iSign -= 2){
- for (int iFace = 0; iFace < 4; iFace++){
- nodalBasis::closure closure_face;
- getFaceClosureTet(iFace, iSign, iRotate, closure_face, order);
- closure.push_back(closure_face);
+ void generateFaceClosureTet(nodalBasis::clCont &closure, int order)
+ {
+ closure.clear();
+ for (int iRotate = 0; iRotate < 3; iRotate++){
+ for (int iSign = 1; iSign >= -1; iSign -= 2){
+ for (int iFace = 0; iFace < 4; iFace++){
+ nodalBasis::closure closure_face;
+ getFaceClosureTet(iFace, iSign, iRotate, closure_face, order);
+ closure.push_back(closure_face);
+ }
}
}
}
-}
-
-
-static void generateFaceClosureTetFull(nodalBasis::clCont &closureFull,
- std::vector<int> &closureRef,
- int order, bool serendip)
-{
- closureFull.clear();
- //input :
- static const short int faces[4][3] = {{0,1,2}, {0,3,1}, {0,2,3}, {3,2,1}};
- static const int edges[] = {0, 1, 1, 2, 2, 0, 3, 0, 3, 2, 3, 1, -1};
- static const int faceOrientation[6] = {0,1,2,5,3,4};
- closureFull.resize(24);
- closureRef.resize(24);
- for (int i = 0; i < 24; i ++)
- closureRef[i] = 0;
- if (order == 0) {
- for (int i = 0; i < 24; i ++) {
- closureFull[i].push_back(0);
- }
- return;
- }
- //Mapping for the p1 nodes
- nodalBasis::clCont closure;
- generateFaceClosureTet(closure, 1);
- 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 (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)
- clFull[j] = (6 - clFull[(j + 1) % 4] - clFull[(j + 2) % 4] - clFull[(j + 3) % 4]);
- }
- int nodes2Faces[4][4][4];
- for (int i = 0; i < 4; i++) {
- for (int iRotate = 0; iRotate < 3; iRotate ++) {
- short int n0 = faces[i][(3 - iRotate) % 3];
- short int n1 = faces[i][(4 - iRotate) % 3];
- short int n2 = faces[i][(5 - iRotate) % 3];
- nodes2Faces[n0][n1][n2] = i * 6 + iRotate;
- nodes2Faces[n0][n2][n1] = i * 6 + iRotate + 3;
+ void generateFaceClosureTetFull(nodalBasis::clCont &closureFull,
+ std::vector<int> &closureRef,
+ int order, bool serendip)
+ {
+ closureFull.clear();
+ //input :
+ static const short int faces[4][3] = {{0,1,2}, {0,3,1}, {0,2,3}, {3,2,1}};
+ static const int edges[] = {0, 1, 1, 2, 2, 0, 3, 0, 3, 2, 3, 1, -1};
+ static const int faceOrientation[6] = {0,1,2,5,3,4};
+ closureFull.resize(24);
+ closureRef.resize(24);
+ for (int i = 0; i < 24; i ++)
+ closureRef[i] = 0;
+ if (order == 0) {
+ for (int i = 0; i < 24; i ++) {
+ closureFull[i].push_back(0);
+ }
+ return;
}
- }
- nodalBasis::clCont closureTriangles;
- std::vector<int> closureTrianglesRef;
- if (order >= 3)
- generate2dEdgeClosureFull(closureTriangles, closureTrianglesRef, order - 3, 3, false);
- addEdgeNodes(closureFull, edges, order);
- for (unsigned int iClosure = 0; iClosure < closureFull.size(); iClosure++) {
- //faces
- std::vector<int> &cl = closureFull[iClosure];
- if (order >= 3) {
- for (int iFace = 0; iFace < 4; iFace ++) {
- int n0 = cl[faces[iFace][0]];
- int n1 = cl[faces[iFace][1]];
- int n2 = cl[faces[iFace][2]];
- short int id = nodes2Faces[n0][n1][n2];
- short int iTriClosure = faceOrientation [id % 6];
- 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]);
+ //Mapping for the p1 nodes
+ nodalBasis::clCont closure;
+ generateFaceClosureTet(closure, 1);
+ 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 (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)
+ clFull[j] = (6 - clFull[(j + 1) % 4] - clFull[(j + 2) % 4] - clFull[(j + 3) % 4]);
+ }
+ int nodes2Faces[4][4][4];
+ for (int i = 0; i < 4; i++) {
+ for (int iRotate = 0; iRotate < 3; iRotate ++) {
+ short int n0 = faces[i][(3 - iRotate) % 3];
+ short int n1 = faces[i][(4 - iRotate) % 3];
+ short int n2 = faces[i][(5 - iRotate) % 3];
+ nodes2Faces[n0][n1][n2] = i * 6 + iRotate;
+ nodes2Faces[n0][n2][n1] = i * 6 + iRotate + 3;
+ }
+ }
+ nodalBasis::clCont closureTriangles;
+ std::vector<int> closureTrianglesRef;
+ if (order >= 3)
+ generate2dEdgeClosureFull(closureTriangles, closureTrianglesRef, order - 3, 3, false);
+ addEdgeNodes(closureFull, edges, order);
+ for (unsigned int iClosure = 0; iClosure < closureFull.size(); iClosure++) {
+ //faces
+ std::vector<int> &cl = closureFull[iClosure];
+ if (order >= 3) {
+ for (int iFace = 0; iFace < 4; iFace ++) {
+ int n0 = cl[faces[iFace][0]];
+ int n1 = cl[faces[iFace][1]];
+ int n2 = cl[faces[iFace][2]];
+ short int id = nodes2Faces[n0][n1][n2];
+ short int iTriClosure = faceOrientation [id % 6];
+ 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]);
+ }
}
}
}
- }
- if (order >= 4 && !serendip) {
- nodalBasis::clCont insideClosure;
- std::vector<int> fakeClosureRef;
- generateFaceClosureTetFull(insideClosure, fakeClosureRef, order - 4, false);
- 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);
+ if (order >= 4 && !serendip) {
+ nodalBasis::clCont insideClosure;
+ std::vector<int> fakeClosureRef;
+ generateFaceClosureTetFull(insideClosure, fakeClosureRef, order - 4, false);
+ 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);
+ }
}
}
-}
-
-
-void rotateHex(int iFace, int iRot, int iSign, double uI, double vI,
- double &uO, double &vO, double &wO)
-{
- if (iSign<0){
- double tmp=uI;
- uI=vI;
- vI=tmp;
- }
- for (int i=0; i < iRot; i++){
- double tmp=uI;
- uI=-vI;
- vI=tmp;
- }
- switch (iFace) {
- case 0: uO = vI; vO = uI; wO=-1; break;
- case 1: uO = uI; vO = -1; wO=vI; break;
- case 2: uO =-1; vO = vI; wO=uI; break;
- case 3: uO = 1; vO = uI; wO=vI; break;
- case 4: uO =-uI; vO = 1; wO=vI; break;
- case 5: uO =uI; vO = vI; wO=1; break;
- }
-}
-
-
-
-void rotateHexFull(int iFace, int iRot, int iSign, double uI, double vI, double wI, double &uO, double &vO, double &wO)
-{
- switch (iFace) {
- case 0: uO = uI; vO = vI; wO = wI; break;
- case 1: uO = wI; vO = uI; wO = vI; break;
- case 2: uO = vI; vO = wI; wO = uI; break;
- case 3: uO = wI; vO = vI; wO =-uI; break;
- case 4: uO = wI; vO =-uI; wO =-vI; break;
- case 5: uO = vI; vO = uI; wO =-wI; break;
- }
- for (int i = 0; i < iRot; i++){
- double tmp = uO;
- uO = -vO;
- vO = tmp;
- }
- if (iSign<0){
- double tmp = uO;
- uO = vO;
- vO = tmp;
- }
-}
-
-
-
-inline static double pow2(double x)
-{
- return x * x;
-}
-
-/*
-static void checkClosure(int tag) {
- printf("TYPE = %i\n", tag);
- const nodalBasis &fs = *nodalBases::find(tag);
- for (int i = 0; i < fs.closures.size(); ++i) {
- const std::vector<int> &clRef = fs.closures[fs.closureRef[i]];
- const std::vector<int> &cl = fs.closures[i];
- const std::vector<int> &clFull = fs.fullClosures[i];
- printf("i = %i\n", i);
- for (int j = 0; j < cl.size(); ++j) {
- printf("%3i ", clFull[clRef[j]]);
- }
- printf("\n");
- for (int j = 0; j < cl.size(); ++j) {
- printf("%3i ", cl[j]);
- }
- printf("\n");
+ /*
+ void checkClosure(int tag) {
+ printf("TYPE = %i\n", tag);
+ const nodalBasis &fs = *nodalBases::find(tag);
+ for (int i = 0; i < fs.closures.size(); ++i) {
+ const std::vector<int> &clRef = fs.closures[fs.closureRef[i]];
+ const std::vector<int> &cl = fs.closures[i];
+ const std::vector<int> &clFull = fs.fullClosures[i];
+ printf("i = %i\n", i);
+ for (int j = 0; j < cl.size(); ++j) {
+ printf("%3i ", clFull[clRef[j]]);
+ }
+ printf("\n");
+ for (int j = 0; j < cl.size(); ++j) {
+ printf("%3i ", cl[j]);
+ }
+ printf("\n");
+ }
}
-}
-*/
-
-
-
-static void generateFaceClosureHex(nodalBasis::clCont &closure, int order,
- bool serendip, const fullMatrix<double> &points)
-{
- closure.clear();
- const nodalBasis &fsFace = *BasisFactory::getNodalBasis(ElementType::getTag(TYPE_QUA, order, serendip));
- for (int iRotate = 0; iRotate < 4; iRotate++){
- for (int iSign = 1; iSign >= -1; iSign -= 2){
- for (int iFace = 0; iFace < 6; iFace++) {
- nodalBasis::closure cl;
- cl.type = fsFace.type;
- cl.resize(fsFace.points.size1());
- for (unsigned int iNode = 0; iNode < cl.size(); ++iNode) {
- double u,v,w;
- rotateHex(iFace, iRotate, iSign, fsFace.points(iNode, 0),
- fsFace.points(iNode, 1), u, v, w);
- cl[iNode] = 0;
- double D = std::numeric_limits<double>::max();
- for (int jNode = 0; jNode < points.size1(); ++jNode) {
- double d = pow2(points(jNode, 0) - u) + pow2(points(jNode, 1) - v) +
- pow2(points(jNode, 2) - w);
- if (d < D) {
- cl[iNode] = jNode;
- D = d;
+ */
+
+ void generateFaceClosureHex(nodalBasis::clCont &closure, int order,
+ bool serendip, const fullMatrix<double> &points)
+ {
+ closure.clear();
+ const nodalBasis &fsFace = *BasisFactory::getNodalBasis(ElementType::getTag(TYPE_QUA, order, serendip));
+ for (int iRotate = 0; iRotate < 4; iRotate++){
+ for (int iSign = 1; iSign >= -1; iSign -= 2){
+ for (int iFace = 0; iFace < 6; iFace++) {
+ nodalBasis::closure cl;
+ cl.type = fsFace.type;
+ cl.resize(fsFace.points.size1());
+ for (unsigned int iNode = 0; iNode < cl.size(); ++iNode) {
+ double u,v,w;
+ rotateHex(iFace, iRotate, iSign, fsFace.points(iNode, 0),
+ fsFace.points(iNode, 1), u, v, w);
+ cl[iNode] = 0;
+ double D = std::numeric_limits<double>::max();
+ for (int jNode = 0; jNode < points.size1(); ++jNode) {
+ double d = pow2(points(jNode, 0) - u) + pow2(points(jNode, 1) - v) +
+ pow2(points(jNode, 2) - w);
+ if (d < D) {
+ cl[iNode] = jNode;
+ D = d;
+ }
}
}
+ closure.push_back(cl);
}
- closure.push_back(cl);
}
}
}
-}
-
-
-static void generateFaceClosureHexFull(nodalBasis::clCont &closure,
- std::vector<int> &closureRef,
- int order, bool serendip,
- const fullMatrix<double> &points)
-{
- closure.clear();
- int clId = 0;
- for (int iRotate = 0; iRotate < 4; iRotate++){
- for (int iSign = 1; iSign >= -1; iSign -= 2){
- for (int iFace = 0; iFace < 6; iFace++) {
- nodalBasis::closure cl;
- cl.resize(points.size1());
- for (int iNode = 0; iNode < points.size1(); ++iNode) {
- double u,v,w;
- rotateHexFull(iFace, iRotate, iSign, points(iNode, 0), points(iNode, 1),
- points(iNode, 2), u, v, w);
- int J = 0;
- double D = std::numeric_limits<double>::max();
- for (int jNode = 0; jNode < points.size1(); ++jNode) {
- double d = pow2(points(jNode, 0) - u) + pow2(points(jNode, 1) - v) +
- pow2(points(jNode, 2) - w);
- if (d < D) {
- J = jNode;
- D = d;
+ void generateFaceClosureHexFull(nodalBasis::clCont &closure,
+ std::vector<int> &closureRef,
+ int order, bool serendip,
+ const fullMatrix<double> &points)
+ {
+ closure.clear();
+ int clId = 0;
+ for (int iRotate = 0; iRotate < 4; iRotate++){
+ for (int iSign = 1; iSign >= -1; iSign -= 2){
+ for (int iFace = 0; iFace < 6; iFace++) {
+ nodalBasis::closure cl;
+ cl.resize(points.size1());
+ for (int iNode = 0; iNode < points.size1(); ++iNode) {
+ double u,v,w;
+ rotateHexFull(iFace, iRotate, iSign, points(iNode, 0), points(iNode, 1),
+ points(iNode, 2), u, v, w);
+ int J = 0;
+ double D = std::numeric_limits<double>::max();
+ for (int jNode = 0; jNode < points.size1(); ++jNode) {
+ double d = pow2(points(jNode, 0) - u) + pow2(points(jNode, 1) - v) +
+ pow2(points(jNode, 2) - w);
+ if (d < D) {
+ J = jNode;
+ D = d;
+ }
}
+ cl[J] = iNode;
}
- cl[J] = iNode;
+ closure.push_back(cl);
+ closureRef.push_back(0);
+ clId ++;
}
- closure.push_back(cl);
- closureRef.push_back(0);
- clId ++;
}
}
}
-}
-
-
-static void getFaceClosurePrism(int iFace, int iSign, int iRotate,
- nodalBasis::closure &closure, int order)
-{
- //if (order > 2)
- // Msg::Error("FaceClosure not implemented for prisms of order %d",order);
- bool isTriangle = iFace<2;
- int nNodes = isTriangle ? (order+1)*(order+2)/2 : (order+1)*(order+1);
- closure.clear();
- if (isTriangle && iRotate > 2) return;
- closure.resize(nNodes);
- if (order==0) {
- closure[0] = 0;
- return;
- }
- int order1node[5][4] = {{0, 2, 1, -1}, {3, 4, 5, -1}, {0, 1, 4, 3}, {0, 3, 5, 2},
- {1, 2, 5, 4}};
- int order2node[5][5] = {{7, 9, 6, -1, -1}, {12, 14, 13, -1, -1}, {6, 10, 12, 8, 15},
- {8, 13, 11, 7, 16}, {9, 11, 14, 10, 17}};
- // int order2node[5][4] = {{7, 9, 6, -1}, {12, 14, 13, -1}, {6, 10, 12, 8},
- // {8, 13, 11, 7}, {9, 11, 14, 10}};
- int nVertex = isTriangle ? 3 : 4;
- closure.type = ElementType::getTag(isTriangle ? TYPE_TRI : TYPE_QUA, order);
- for (int i = 0; i < nVertex; ++i){
- int k = (nVertex + (iSign * i) + iRotate) % nVertex; //- iSign * iRotate
- closure[i] = order1node[iFace][k];
- }
- if (order==2) {
+ void getFaceClosurePrism(int iFace, int iSign, int iRotate,
+ nodalBasis::closure &closure, int order)
+ {
+ //if (order > 2)
+ // Msg::Error("FaceClosure not implemented for prisms of order %d",order);
+ bool isTriangle = iFace<2;
+ int nNodes = isTriangle ? (order+1)*(order+2)/2 : (order+1)*(order+1);
+ closure.clear();
+ if (isTriangle && iRotate > 2) return;
+ closure.resize(nNodes);
+ if (order==0) {
+ closure[0] = 0;
+ return;
+ }
+ int order1node[5][4] = {{0, 2, 1, -1}, {3, 4, 5, -1}, {0, 1, 4, 3}, {0, 3, 5, 2},
+ {1, 2, 5, 4}};
+ int order2node[5][5] = {{7, 9, 6, -1, -1}, {12, 14, 13, -1, -1}, {6, 10, 12, 8, 15},
+ {8, 13, 11, 7, 16}, {9, 11, 14, 10, 17}};
+ // int order2node[5][4] = {{7, 9, 6, -1}, {12, 14, 13, -1}, {6, 10, 12, 8},
+ // {8, 13, 11, 7}, {9, 11, 14, 10}};
+ int nVertex = isTriangle ? 3 : 4;
+ closure.type = ElementType::getTag(isTriangle ? TYPE_TRI : TYPE_QUA, order);
for (int i = 0; i < nVertex; ++i){
- int k = (nVertex + (iSign==-1?-1:0) + (iSign * i) + iRotate) % nVertex;
- //- iSign * iRotate
- closure[nVertex+i] = order2node[iFace][k];
+ int k = (nVertex + (iSign * i) + iRotate) % nVertex; //- iSign * iRotate
+ closure[i] = order1node[iFace][k];
}
- if (!isTriangle)
- closure[nNodes-1] = order2node[iFace][4]; // center
- }
-}
-
-
-
-static void generateFaceClosurePrism(nodalBasis::clCont &closure, int order)
-{
- if (order > 2)
- Msg::Error("FaceClosure not implemented for prisms of order %d",order);
- closure.clear();
- for (int iRotate = 0; iRotate < 4; iRotate++){
- for (int iSign = 1; iSign >= -1; iSign -= 2){
- for (int iFace = 0; iFace < 5; iFace++){
- nodalBasis::closure closure_face;
- getFaceClosurePrism(iFace, iSign, iRotate, closure_face, order);
- closure.push_back(closure_face);
+ if (order==2) {
+ for (int i = 0; i < nVertex; ++i){
+ int k = (nVertex + (iSign==-1?-1:0) + (iSign * i) + iRotate) % nVertex;
+ //- iSign * iRotate
+ closure[nVertex+i] = order2node[iFace][k];
}
+ if (!isTriangle)
+ closure[nNodes-1] = order2node[iFace][4]; // center
}
}
-}
-
-
-static void generateFaceClosurePrismFull(nodalBasis::clCont &closureFull,
- std::vector<int> &closureRef, int order)
-{
- nodalBasis::clCont closure;
- closureFull.clear();
- closureFull.resize(40);
- closureRef.resize(40);
- generateFaceClosurePrism(closure, 1);
- int ref3 = -1, ref4a = -1, ref4b = -1;
- for (unsigned int i = 0; i < closure.size(); i++) {
- std::vector<int> &clFull = closureFull[i];
- std::vector<int> &cl = closure[i];
- if (cl.size() == 0)
- continue;
- clFull.resize(6, -1);
- int &ref = cl.size() == 3 ? ref3 : (cl[0] / 3 + cl[1] / 3) % 2 ? ref4b : ref4a;
- if (ref == -1)
- ref = i;
- closureRef[i] = ref;
- 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) {
- int k = ((j / 3) + 1) % 2 * 3;
- int sum = (clFull[k + (j + 1) % 3] + clFull[k + (j + 2) % 3]);
- clFull[j] = ((sum / 6 + 1) % 2) * 3 + (12 - sum) % 3;
+ void generateFaceClosurePrism(nodalBasis::clCont &closure, int order)
+ {
+ if (order > 2)
+ Msg::Error("FaceClosure not implemented for prisms of order %d",order);
+ closure.clear();
+ for (int iRotate = 0; iRotate < 4; iRotate++){
+ for (int iSign = 1; iSign >= -1; iSign -= 2){
+ for (int iFace = 0; iFace < 5; iFace++){
+ nodalBasis::closure closure_face;
+ getFaceClosurePrism(iFace, iSign, iRotate, closure_face, order);
+ closure.push_back(closure_face);
+ }
}
}
}
- 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}};
- if ( order == 2 ) {
- int nextFaceNode = 15;
- int numFaces = 5;
- int numFaceNodes = 4;
- std::map<int,int> nodeSum2Face;
- for (int iFace = 0; iFace < numFaces ; iFace ++) {
- int nodeSum = 0;
- for (int iNode = 0; iNode < numFaceNodes; iNode++ ) {
- nodeSum += faces[iFace][iNode];
+ void generateFaceClosurePrismFull(nodalBasis::clCont &closureFull,
+ std::vector<int> &closureRef, int order)
+ {
+ nodalBasis::clCont closure;
+ closureFull.clear();
+ closureFull.resize(40);
+ closureRef.resize(40);
+ generateFaceClosurePrism(closure, 1);
+ int ref3 = -1, ref4a = -1, ref4b = -1;
+ for (unsigned int i = 0; i < closure.size(); i++) {
+ std::vector<int> &clFull = closureFull[i];
+ std::vector<int> &cl = closure[i];
+ if (cl.size() == 0)
+ continue;
+ clFull.resize(6, -1);
+ int &ref = cl.size() == 3 ? ref3 : (cl[0] / 3 + cl[1] / 3) % 2 ? ref4b : ref4a;
+ if (ref == -1)
+ ref = i;
+ closureRef[i] = ref;
+ 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) {
+ int k = ((j / 3) + 1) % 2 * 3;
+ int sum = (clFull[k + (j + 1) % 3] + clFull[k + (j + 2) % 3]);
+ clFull[j] = ((sum / 6 + 1) % 2) * 3 + (12 - sum) % 3;
+ }
}
- nodeSum2Face[nodeSum] = iFace;
}
- for (unsigned int i = 0; i < closureFull.size(); i++ ) {
- if (closureFull[i].empty())
- continue;
- for (int iFace = 0; iFace < numFaces; iFace++ ) {
+ 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}};
+
+ if ( order == 2 ) {
+ int nextFaceNode = 15;
+ int numFaces = 5;
+ int numFaceNodes = 4;
+ std::map<int,int> nodeSum2Face;
+ 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] ];
- std::map<int,int>::iterator it = nodeSum2Face.find(nodeSum);
- if (it == nodeSum2Face.end() )
- Msg::Error("Prism face not found");
- int mappedFaceId = it->second;
- if ( mappedFaceId > 1) {
- closureFull[i].push_back(nextFaceNode + mappedFaceId - 2);
+ for (int iNode = 0; iNode < numFaceNodes; iNode++ ) {
+ nodeSum += faces[iFace][iNode];
}
+ nodeSum2Face[nodeSum] = iFace;
}
+ for (unsigned int i = 0; i < closureFull.size(); i++ ) {
+ if (closureFull[i].empty())
+ continue;
+ 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] ];
+ std::map<int,int>::iterator it = nodeSum2Face.find(nodeSum);
+ if (it == nodeSum2Face.end() )
+ Msg::Error("Prism face not found");
+ int mappedFaceId = it->second;
+ if ( mappedFaceId > 1) {
+ closureFull[i].push_back(nextFaceNode + mappedFaceId - 2);
+ }
+ }
+ }
+ } else {
+ Msg::Error("FaceClosureFull not implemented for prisms of order %d",order);
}
- } else {
- Msg::Error("FaceClosureFull not implemented for prisms of order %d",order);
- }
-}
-
-
-
-static void generate2dEdgeClosure(nodalBasis::clCont &closure, int order, int nNod = 3)
-{
- closure.clear();
- closure.resize(2*nNod);
- for (int j = 0; j < nNod ; j++){
- closure[j].push_back(j);
- closure[j].push_back((j+1)%nNod);
- closure[nNod+j].push_back((j+1)%nNod);
- closure[nNod+j].push_back(j);
- for (int i=0; i < order-1; i++){
- closure[j].push_back( nNod + (order-1)*j + i );
- closure[nNod+j].push_back(nNod + (order-1)*(j+1) -i -1);
- }
- closure[j].type = closure[nNod+j].type = ElementType::getTag(TYPE_LIN, order);
}
-}
-
+ void generate2dEdgeClosure(nodalBasis::clCont &closure, int order, int nNod = 3)
+ {
+ closure.clear();
+ closure.resize(2*nNod);
+ for (int j = 0; j < nNod ; j++){
+ closure[j].push_back(j);
+ closure[j].push_back((j+1)%nNod);
+ closure[nNod+j].push_back((j+1)%nNod);
+ closure[nNod+j].push_back(j);
+ for (int i=0; i < order-1; i++){
+ closure[j].push_back( nNod + (order-1)*j + i );
+ closure[nNod+j].push_back(nNod + (order-1)*(j+1) -i -1);
+ }
+ closure[j].type = closure[nNod+j].type = ElementType::getTag(TYPE_LIN, order);
+ }
+ }
-static void generateClosureOrder0(nodalBasis::clCont &closure, int nb)
-{
- closure.clear();
- closure.resize(nb);
- for (int i=0; i < nb; i++) {
- closure[i].push_back(0);
- closure[i].type = MSH_PNT;
+ void generateClosureOrder0(nodalBasis::clCont &closure, int nb)
+ {
+ closure.clear();
+ closure.resize(nb);
+ for (int i=0; i < nb; i++) {
+ closure[i].push_back(0);
+ closure[i].type = MSH_PNT;
+ }
}
}
-
-
nodalBasis::nodalBasis(int tag)
{
+ using namespace ClosureGen;
type = tag;
parentType = ElementType::ParentTypeFromTag(tag);
order = ElementType::OrderFromTag(tag);
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 3f985904f18b423f3c382839b1fd4c2d9fd7019f..01388992be9c09aa132c59a2658d0bcc1953b74a 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -17,60 +17,61 @@
#include "pointsGenerators.h"
#include <limits>
-/*
-static void printClosure(polynomialBasis::clCont &fullClosure,
- std::vector<int> &closureRef,
- polynomialBasis::clCont &closures)
-{
- for (unsigned int i = 0; i < closures.size(); i++) {
- printf("%3i (%2i): ", i, closureRef[i]);
- if(closureRef[i]==-1){
- printf("\n");
- continue;
+namespace {
+ fullMatrix<double> generateLagrangeMonomialCoefficients
+ (const fullMatrix<double>& monomial, const fullMatrix<double>& point)
+ {
+ if(monomial.size1() != point.size1() || monomial.size2() != point.size2()){
+ Msg::Fatal("Wrong sizes for Lagrange coefficients generation %d %d -- %d %d",
+ monomial.size1(),point.size1(),
+ monomial.size2(),point.size2() );
+ return fullMatrix<double>(1, 1);
}
- for (unsigned int j = 0; j < fullClosure[i].size(); j++) {
- printf("%2i ", fullClosure[i][j]);
- }
- printf (" -- ");
- 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]);
+
+ int ndofs = monomial.size1();
+ int dim = monomial.size2();
+
+ fullMatrix<double> Vandermonde(ndofs, ndofs);
+ for (int i = 0; i < ndofs; i++) {
+ for (int j = 0; j < ndofs; j++) {
+ double dd = 1.;
+ for (int k = 0; k < dim; k++) dd *= pow_int(point(j, k), monomial(i, k));
+ Vandermonde(i, j) = dd;
+ }
}
- printf("\n");
- }
-}
-*/
+ fullMatrix<double> coefficient(ndofs, ndofs);
+ Vandermonde.invert(coefficient);
-static fullMatrix<double> generateLagrangeMonomialCoefficients
- (const fullMatrix<double>& monomial, const fullMatrix<double>& point)
-{
- if(monomial.size1() != point.size1() || monomial.size2() != point.size2()){
- Msg::Fatal("Wrong sizes for Lagrange coefficients generation %d %d -- %d %d",
- monomial.size1(),point.size1(),
- monomial.size2(),point.size2() );
- return fullMatrix<double>(1, 1);
+ return coefficient;
}
- int ndofs = monomial.size1();
- int dim = monomial.size2();
-
- fullMatrix<double> Vandermonde(ndofs, ndofs);
- for (int i = 0; i < ndofs; i++) {
- for (int j = 0; j < ndofs; j++) {
- double dd = 1.;
- for (int k = 0; k < dim; k++) dd *= pow_int(point(j, k), monomial(i, k));
- Vandermonde(i, j) = dd;
+ /*
+ void printClosure(polynomialBasis::clCont &fullClosure,
+ std::vector<int> &closureRef,
+ polynomialBasis::clCont &closures)
+ {
+ for (unsigned int i = 0; i < closures.size(); i++) {
+ printf("%3i (%2i): ", i, closureRef[i]);
+ if(closureRef[i]==-1){
+ printf("\n");
+ continue;
+ }
+ for (unsigned int j = 0; j < fullClosure[i].size(); j++) {
+ printf("%2i ", fullClosure[i][j]);
+ }
+ printf (" -- ");
+ 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("\n");
}
}
-
- fullMatrix<double> coefficient(ndofs, ndofs);
- Vandermonde.invert(coefficient);
-
- return coefficient;
+ */
}
polynomialBasis::polynomialBasis(int tag) : nodalBasis(tag)
@@ -105,15 +106,10 @@ polynomialBasis::polynomialBasis(int tag) : nodalBasis(tag)
coefficients = generateLagrangeMonomialCoefficients(monomials, points);
}
-
-
polynomialBasis::~polynomialBasis()
{
}
-
-
-
void polynomialBasis::evaluateMonomials(double u, double v, double w, double p[]) const
{
for (int j = 0; j < monomials.size1(); ++j) {
@@ -128,8 +124,6 @@ void polynomialBasis::evaluateMonomials(double u, double v, double w, double p[]
}
}
-
-
void polynomialBasis::f(double u, double v, double w, double *sf) const
{
double p[1256];
@@ -142,7 +136,6 @@ void polynomialBasis::f(double u, double v, double w, double *sf) const
}
}
-
void polynomialBasis::f(const fullMatrix<double> &coord, fullMatrix<double> &sf) const
{
double p[1256];
@@ -156,8 +149,6 @@ void polynomialBasis::f(const fullMatrix<double> &coord, fullMatrix<double> &sf)
}
}
-
-
void polynomialBasis::df(const fullMatrix<double> &coord, fullMatrix<double> &dfm) const
{
double dfv[1256][3];
@@ -173,8 +164,6 @@ void polynomialBasis::df(const fullMatrix<double> &coord, fullMatrix<double> &df
}
}
-
-
void polynomialBasis::df(double u, double v, double w, double grads[][3]) const
{
switch (monomials.size2()) {
@@ -234,8 +223,6 @@ void polynomialBasis::df(double u, double v, double w, double grads[][3]) const
}
}
-
-
void polynomialBasis::ddf(double u, double v, double w, double hess[][3][3]) const
{
switch (monomials.size2()) {
@@ -319,8 +306,6 @@ void polynomialBasis::ddf(double u, double v, double w, double hess[][3][3]) con
}
}
-
-
void polynomialBasis::dddf(double u, double v, double w, double third[][3][3][3]) const
{
switch (monomials.size2()) {
diff --git a/Plugin/AnalyseCurvedMesh.cpp b/Plugin/AnalyseCurvedMesh.cpp
index 58f3d04685dd8d4d79fca6f8bbf09acf1719aca8..b73f262befe0e1db51790f0aec82f6891bca8e97 100644
--- a/Plugin/AnalyseCurvedMesh.cpp
+++ b/Plugin/AnalyseCurvedMesh.cpp
@@ -23,6 +23,17 @@
#include "FlGui.h"
#endif
+namespace {
+ double sum(fullVector<double> &v)
+ {
+ double sum = .0;
+ for (int i = 0; i < v.size(); ++i) {
+ sum += v(i);
+ }
+ return sum;
+ }
+}
+
//#define UNDEF_JAC_TAG -999
//#define _ANALYSECURVEDMESH_BLAS_
@@ -75,15 +86,6 @@ std::string GMSH_AnalyseCurvedMeshPlugin::getHelp() const
"Tolerance = R+ , << 1 : tolerance (relatively to J_min and J_max) used during the computation of J_min and J_max";
}
-static double sum(fullVector<double> &v)
-{
- double sum = .0;
- for (int i = 0; i < v.size(); ++i) {
- sum += v(i);
- }
- return sum;
-}
-
// Execution
PView *GMSH_AnalyseCurvedMeshPlugin::execute(PView *v)
{