diff --git a/Geo/MElementCut.h b/Geo/MElementCut.h
index d908a0cd2f4dc906c1c0e8d5fc69ced5300c9705..137abb1615248c3c71adec15d6e2045005e19898 100644
--- a/Geo/MElementCut.h
+++ b/Geo/MElementCut.h
@@ -119,6 +119,11 @@ class MPolyhedron : public MElement {
if(_orig) return _orig->getFunctionSpace(order);
return 0;
}
+ virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const
+ {
+ if(_orig) return _orig->getJacobianFuncSpace(order);
+ return 0;
+ }
virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
{
if(_orig) _orig->getShapeFunctions(u, v, w, s, o);
@@ -241,6 +246,11 @@ class MPolygon : public MElement {
if(_orig) return _orig->getFunctionSpace(order);
return 0;
}
+ virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const
+ {
+ if(_orig) return _orig->getJacobianFuncSpace(order);
+ return 0;
+ }
virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
{
if(_orig) _orig->getShapeFunctions(u, v, w, s, o);
@@ -297,6 +307,11 @@ class MLineChild : public MLine {
if(_orig) return _orig->getFunctionSpace(order);
return 0;
}
+ virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const
+ {
+ if(_orig) return _orig->getJacobianFuncSpace(order);
+ return 0;
+ }
virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
{
if(_orig) _orig->getShapeFunctions(u, v, w, s, o);
diff --git a/Geo/MLine.cpp b/Geo/MLine.cpp
index 36e1bb5cdf2152a688f916603e96c65ad7e36263..5eaed4c83651c2b145dc7fdf314809483c0046de 100644
--- a/Geo/MLine.cpp
+++ b/Geo/MLine.cpp
@@ -29,6 +29,26 @@ const polynomialBasis* MLine::getFunctionSpace(int o) const
return 0;
}
+const JacobianBasis* MLine::getJacobianFuncSpace(int o) const
+{
+ int order = (o == -1) ? getPolynomialOrder() : o;
+
+ switch (order) {
+ case 1: return JacobianBases::find(MSH_LIN_2);
+ case 2: return JacobianBases::find(MSH_LIN_3);
+ case 3: return JacobianBases::find(MSH_LIN_4);
+ case 4: return JacobianBases::find(MSH_LIN_5);
+ case 5: return JacobianBases::find(MSH_LIN_6);
+ case 6: return JacobianBases::find(MSH_LIN_7);
+ case 7: return JacobianBases::find(MSH_LIN_8);
+ case 8: return JacobianBases::find(MSH_LIN_9);
+ case 9: return JacobianBases::find(MSH_LIN_10);
+ case 10: return JacobianBases::find(MSH_LIN_11);
+ default: Msg::Error("Order %d line function space not implemented", order);
+ }
+ return 0;
+}
+
void MLine::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
{
*npts = getNGQLPts(pOrder);
diff --git a/Geo/MLine.h b/Geo/MLine.h
index 567414087f0d1458c72c271c81d4445bf7e8c43f..0f66d608b13fbb8a17b6e048e0aa8bcbb2f2d1f8 100644
--- a/Geo/MLine.h
+++ b/Geo/MLine.h
@@ -70,6 +70,7 @@ class MLine : public MElement {
MVertex *tmp = _v[0]; _v[0] = _v[1]; _v[1] = tmp;
}
virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
+ virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
virtual bool isInside(double u, double v, double w)
{
double tol = _isInsideTolerance;
diff --git a/Geo/MPoint.h b/Geo/MPoint.h
index 81b6f2cb3dd0ee3033b712a06a93ac8e3a9da678..993e7743358b532798fce65d8cb9a8395b2fe63d 100644
--- a/Geo/MPoint.h
+++ b/Geo/MPoint.h
@@ -52,6 +52,7 @@ class MPoint : public MElement {
}
virtual const polynomialBasis* getFunctionSpace(int o) const { return polynomialBases::find(MSH_PNT); }
+ virtual const JacobianBasis* getJacobianFuncSpace(int o) const { return JacobianBases::find(MSH_PNT); }
virtual bool isInside(double u, double v, double w)
{
return true;
diff --git a/Geo/MPrism.cpp b/Geo/MPrism.cpp
index d6e99b70637c3b896acfbb0e3bb41378a9721bbf..8333b2e576a4da9f860f6838d74a51012b573b86 100644
--- a/Geo/MPrism.cpp
+++ b/Geo/MPrism.cpp
@@ -53,6 +53,29 @@ const polynomialBasis* MPrism::getFunctionSpace(int o) const
return 0;
}
+const JacobianBasis* MPrism::getJacobianFuncSpace(int o) const
+{
+ int order = (o == -1) ? getPolynomialOrder() : o;
+
+ int nv = getNumVolumeVertices();
+
+ if ((nv == 0) && (o == -1)) {
+ switch (order) {
+ case 1: return JacobianBases::find(MSH_PRI_6);
+ case 2: return JacobianBases::find(MSH_PRI_18);
+ default: Msg::Error("Order %d prism function space not implemented", order);
+ }
+ }
+ else {
+ switch (order) {
+ case 1: return JacobianBases::find(MSH_PRI_6);
+ case 2: return JacobianBases::find(MSH_PRI_18);
+ default: Msg::Error("Order %d prism function space not implemented", order);
+ }
+ }
+ return 0;
+}
+
double MPrism::getInnerRadius()
{
double dist[3], k = 0.;
diff --git a/Geo/MPrism.h b/Geo/MPrism.h
index 9336a490ffd0ed093dfdbc2dbf0cbe5e5f730a99..4703450d5e05feca50b04c0a96c7eda559069f4e 100644
--- a/Geo/MPrism.h
+++ b/Geo/MPrism.h
@@ -133,6 +133,7 @@ class MPrism : public MElement {
tmp = _v[3]; _v[3] = _v[4]; _v[4] = tmp;
}
virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
+ virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
virtual int getVolumeSign();
/* virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
{
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index 92eacb3a1bfd9464fad8520f21af9b0c420a51c9..bc3b8fd609304a72834f1062fa56df12911ee65a 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -52,6 +52,42 @@ const polynomialBasis* MQuadrangle::getFunctionSpace(int o) const
return 0;
}
+const JacobianBasis* MQuadrangle::getJacobianFuncSpace(int o) const
+{
+ int order = (o == -1) ? getPolynomialOrder() : o;
+
+ int nf = getNumFaceVertices();
+
+ if ((nf == 0) && (o == -1)) {
+ switch (order) {
+ case 1: return JacobianBases::find(MSH_QUA_4);
+ case 2: return JacobianBases::find(MSH_QUA_8);
+ case 3: return JacobianBases::find(MSH_QUA_12);
+ case 4: return JacobianBases::find(MSH_QUA_16I);
+ case 5: return JacobianBases::find(MSH_QUA_20);
+ case 6: return JacobianBases::find(MSH_QUA_24);
+ case 7: return JacobianBases::find(MSH_QUA_28);
+ case 8: return JacobianBases::find(MSH_QUA_32);
+ case 9: return JacobianBases::find(MSH_QUA_36I);
+ case 10: return JacobianBases::find(MSH_QUA_40);
+ }
+ }
+ switch (order) {
+ case 1: return JacobianBases::find(MSH_QUA_4);
+ case 2: return JacobianBases::find(MSH_QUA_9);
+ case 3: return JacobianBases::find(MSH_QUA_16);
+ case 4: return JacobianBases::find(MSH_QUA_25);
+ case 5: return JacobianBases::find(MSH_QUA_36);
+ case 6: return JacobianBases::find(MSH_QUA_49);
+ case 7: return JacobianBases::find(MSH_QUA_64);
+ case 8: return JacobianBases::find(MSH_QUA_81);
+ case 9: return JacobianBases::find(MSH_QUA_100);
+ case 10: return JacobianBases::find(MSH_QUA_121);
+ default: Msg::Error("Order %d triangle function space not implemented", order);
+ }
+ return 0;
+}
+
int MQuadrangleN::getNumEdgesRep(){ return 4 * CTX::instance()->mesh.numSubEdges; }
int MQuadrangle8::getNumEdgesRep(){ return 4 * CTX::instance()->mesh.numSubEdges; }
int MQuadrangle9::getNumEdgesRep(){ return 4 * CTX::instance()->mesh.numSubEdges; }
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index ea123cd7f751f51919bf7c87e03fc61bbd6dfc00..4684e2b0c7a1b9012d1596517de83fc1026af2ab 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -119,6 +119,7 @@ class MQuadrangle : public MElement {
virtual const char *getStringForDIFF() const { return "ElmB4n2D"; }
virtual const char *getStringForINP() const { return "C2D4"; }
virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
+ virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
virtual void revert()
{
MVertex *tmp = _v[1]; _v[1] = _v[3]; _v[3] = tmp;
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index 4ff2f8b933c1b9bde6909c1168df5658ae96d1c5..2863cbceac9f345d49375f8b63426953e857a522 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -1,1679 +1,1391 @@
-// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to <gmsh@geuz.org>.
-//
-// Contributor(s):
-// Koen Hillewaert
-//
-
-#include "GmshDefines.h"
-#include "GmshMessage.h"
-#include "JacobianBasis.h"
-
-
-
-#include "polynomialBasis.h"
-
-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;
- }
- exposants.print("1DExp");
- 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);
- }
- }
- }
- exposants.print("ExpTriangle");
- return exposants;
-}
-//
-//// generate the exterior hull of the Exposants triangle for Serendipity element
-//
-//static fullMatrix<double> generateExposantsSerendipityTriangle(int order)
-//{
-// fullMatrix<double> monomials(3 * order, 2);
-//
-// monomials(0, 0) = 0;
-// monomials(0, 1) = 0;
-//
-// int index = 1;
-// for (int i = 1; i <= order; i++) {
-// if (i == order) {
-// for (int j = 0; j <= i; j++) {
-// monomials(index, 0) = i - j;
-// monomials(index, 1) = j;
-// index++;
-// }
-// }
-// else {
-// monomials(index, 0) = i;
-// monomials(index, 1) = 0;
-// index++;
-// monomials(index, 0) = 0;
-// monomials(index, 1) = i;
-// index++;
-// }
-// }
-// return monomials;
-//}
-
-// generate all monomials xi^m * eta^n with n and m <= order
-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);
- }
- }
- }
- exposants.print("ExpQuad");
-
- return exposants;
-}
-///*
-//00 10 20 30 40 ⋯
-//01 11 21 31 41 ⋯
-//02 12
-//03 13
-//04 14
-//â‹® â‹®
-//*/
-//static fullMatrix<double> generateExposantsQuadSerendip(int order)
-//{
-// fullMatrix<double> monomials( (order)*4, 2);
-// monomials(0,0)=0.;
-// monomials(0,1)=0.;
-// monomials(1,0)=1.;
-// monomials(1,1)=0.;
-// monomials(2,0)=0.;
-// monomials(2,1)=1.;
-// monomials(3,0)=1.;
-// monomials(3,1)=1.;
-// int index = 4;
-// for (int p = 2; p <= order; p++) {
-// monomials(index, 0) = p;
-// monomials(index, 1) = 0;
-// index++;
-// monomials(index, 0) = 0;
-// monomials(index, 1) = p;
-// index++;
-// monomials(index, 0) = p;
-// monomials(index, 1) = 1;
-// index++;
-// monomials(index, 0) = 1;
-// monomials(index, 1) = p;
-// index++;
-// }
-// return monomials;
-//}
-//
-///*static fullMatrix<double> generateExposantsQuadSerendip(int order)
-//{
-//
-// fullMatrix<double> monomials( order*4, 2);
-// int index = 0;
-// for (int p = 0; p < order; p++) {
-// monomials(p, 0) = p;
-// monomials(p, 1) = 0;
-//
-// monomials(p+order, 0) = order;
-// monomials(p+order, 1) = p;
-//
-// monomials(p+3*order, 0) = order-p;
-// monomials(p+3*order, 1) = order;
-//
-// monomials(p+2*order, 0) = 0;
-// monomials(p+2*order, 1) = order-p;
-// }
-// monomials.print();
-// return monomials;
-//}*/
-//
-//// generate the monomials subspace of all monomials of order exactly == p
-//
-//static fullMatrix<double> generateMonomialSubspace(int dim, int p)
-//{
-// fullMatrix<double> monomials;
-//
-// switch (dim) {
-// case 1:
-// monomials = fullMatrix<double>(1, 1);
-// monomials(0, 0) = p;
-// break;
-// case 2:
-// monomials = fullMatrix<double>(p + 1, 2);
-// for (int k = 0; k <= p; k++) {
-// monomials(k, 0) = p - k;
-// monomials(k, 1) = k;
-// }
-// break;
-// case 3:
-// monomials = fullMatrix<double>((p + 1) * (p + 2) / 2, 3);
-// int index = 0;
-// for (int i = 0; i <= p; i++) {
-// for (int k = 0; k <= p - i; k++) {
-// monomials(index, 0) = p - i - k;
-// monomials(index, 1) = k;
-// monomials(index, 2) = i;
-// index++;
-// }
-// }
-// break;
-// }
-// return monomials;
-//}
-//
-//// generate external hull of the Exposants tetrahedron
-//
-//static fullMatrix<double> generateExposantsSerendipityTetrahedron(int order)
-//{
-// int nbMonomials = 4 + 6 * std::max(0, order - 1) +
-// 4 * std::max(0, (order - 2) * (order - 1) / 2);
-// fullMatrix<double> monomials(nbMonomials, 3);
-//
-// monomials.setAll(0);
-// int index = 1;
-// for (int p = 1; p < order; p++) {
-// for (int i = 0; i < 3; i++) {
-// int j = (i + 1) % 3;
-// int k = (i + 2) % 3;
-// for (int ii = 0; ii < p; ii++) {
-// monomials(index, i) = p - ii;
-// monomials(index, j) = ii;
-// monomials(index, k) = 0;
-// index++;
-// }
-// }
-// }
-// fullMatrix<double> monomialsMaxOrder = generateMonomialSubspace(3, order);
-// int nbMaxOrder = monomialsMaxOrder.size1();
-// monomials.copy(monomialsMaxOrder, 0, nbMaxOrder, 0, 3, index, 0);
-// return monomials;
-//}
-//
-static int nbdoftriangle(int order) { return (order + 1) * (order + 2) / 2; }
-//static int nbdoftriangleserendip(int order) { return 3 * order; }
-
-//KH : caveat : node coordinates are not yet coherent with node numbering associated
-// to numbering of principal vertices of face !!!!
-
-// uv surface - orientation v0-v2-v1
-static void nodepositionface0(int order, double *u, double *v, double *w)
-{
- 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;
- }
-}
-
-// uw surface - orientation v0-v1-v3
-static void nodepositionface1(int order, double *u, double *v, double *w)
-{
- 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;
- }
-}
-
-// vw surface - orientation v0-v3-v2
-static void nodepositionface2(int order, double *u, double *v, double *w)
-{
- 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;
- }
-}
-
-// uvw surface - orientation v3-v1-v2
-static void nodepositionface3(int order, double *u, double *v, double *w)
-{
- 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;
- }
-}
-
-// generate Exposants tetrahedron
-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);
- }
- }
- }
- }
- }
- exposants.print("ExpTet");
-
- return exposants;
-}
-
-// generate Exposants prism
-
-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*/
- {
- int i, j;
- for (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 (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 ++;
- }
- }
- }
-
-// exposants.print("Pri ipts");
- exposants.print("ExpPrism");
-
- return exposants;
-
-}
-
-static fullMatrix<double> generate1DPoints(int order)
-{
- fullMatrix<double> line(order + 1, 1);
- line(0,0) = 0;
- if (order > 0) {
- line(0, 0) = -1.;
- line(1, 0) = 1.;
- double dd = 2. / order;
- for (int i = 2; i < order + 1; i++) line(i, 0) = -1. + dd * (i - 1);
- }
- line.print("Line");
- 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> gmshGeneratePointsPrism(int order, bool serendip)
-{
- 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> point(nbPoints, 3);
-
- int index = 0;
- fullMatrix<double> triPoints = gmshGeneratePointsTriangle(order,false);
- fullMatrix<double> linePoints = generate1DPoints(order);
-
- if (order == 2)
- for (int i =0; i<18; i++)
- for (int j=0; j<3;j++)
- point(i,j) = prism18Pts[i][j];
- else
- 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 ++;
- }
- }
-
-// point.print("Pri ipts");
-
- return point;
-}
-
-static fullMatrix<double> gmshGeneratePointsQuad(int order, bool serendip)
-{
- int nbPoints = serendip ? order*4 : (order+1)*(order+1);
- fullMatrix<double> point(nbPoints, 2);
-
- if (order > 0) {
- point(0, 0) = -1;
- point(0, 1) = -1;
- point(1, 0) = 1;
- point(1, 1) = -1;
- point(2, 0) = 1;
- point(2, 1) = 1;
- point(3, 0) = -1;
- 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 = gmshGeneratePointsQuad(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;
-}
-
-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;
- // attention, c*n est toujours multiple de i mais
- // n/i n'est pas entier ! Vérifier ordre des opérations.
- return c;
-}
-
-static fullMatrix<double> generateBezierFFCoefficientsSimplex
- (const fullMatrix<double>& monomial, const fullMatrix<double>& point, int order)
-{ // A simplex is a generalization of the notion of a triangle to arbitrary dimension
-
- if(monomial.size1() != point.size1() || monomial.size2() != point.size2()){
- Msg::Fatal("Wrong sizes for Bezier coefficients generation %d %d -- %d %d",
- monomial.size1(),point.size1(),
- monomial.size2(),point.size2());
- return fullMatrix<double>(1, 1);
- }
-
- 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.;
- double pointCompl = 1.;
- int monomialCompl = order;
- for (int k = 0; k < dim; k++) {
- dd *= nChoosek(monomialCompl, (int) monomial(i, k)) * pow(point(j, k), monomial(i, k));
- pointCompl -= point(j, k);
- monomialCompl -= (int) monomial(i, k);
- }
- Vandermonde(j, i) = dd * pow(pointCompl, monomialCompl);
- }
- }
-
- fullMatrix<double> coefficient(ndofs, ndofs);
- Vandermonde.invert(coefficient);
- return coefficient;
-}
-//
-//static void getFaceClosure(int iFace, int iSign, int iRotate, std::vector<int> &closure,
-// int order)
-//{
-//
-// closure.clear();
-// closure.resize((order + 1) * (order + 2) / 2);
-// 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;
-// }
-// }
-// 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++;
-// }
-// }
-// break;
-// }
-//
-//}
-//
-//static void generate3dFaceClosure(polynomialBasis::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++){
-// std::vector<int> closure_face;
-// getFaceClosure(iFace, iSign, iRotate, closure_face, order);
-// closure.push_back(closure_face);
-// }
-// }
-// }
-//}
-//
-//static void getFaceClosurePrism(int iFace, int iSign, int iRotate, std::vector<int> &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();
-// 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;
-// for (int i = 0; i < nVertex; ++i){
-// int k = (nVertex + (iSign * i) + iRotate) % nVertex; //- iSign * iRotate
-// closure[i] = order1node[iFace][k];
-// }
-// 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 generate3dFaceClosurePrism(polynomialBasis::clCont &closure, int 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++){
-// std::vector<int> closure_face;
-// getFaceClosurePrism(iFace, iSign, iRotate, closure_face, order);
-// closure.push_back(closure_face);
-// }
-// }
-// }
-//}
-//
-//
-//static void generate2dEdgeClosure(polynomialBasis::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);
-// }
-// }
-//}
-//
-//static void generate1dVertexClosure(polynomialBasis::clCont &closure)
-//{
-//
-// closure.clear();
-// closure.resize(2);
-// closure[0].push_back(0);
-// closure[1].push_back(1);
-//
-//}
-//
-// F.coefficients = generateLagrangeMonomialCoefficients(F.monomials, F.points);
-//// printf("Case: %d coeffs:\n",tag);
-//// for (int i = 0; i<F.coefficients.size1(); i++) {
-//// for (int j = 0; j<F.coefficients.size2(); j++) {
-//// printf("%4.1f ",F.coefficients(i,j));
-//// }
-//// printf("\n");
-//// }
-//
-// fs.insert(std::make_pair(tag, F));
-// return &fs[tag];
-//}
-
-
-
-// A faire :
-// - getjacobianFunctionSpace for point, line, quad, tet, prism, hexa
-// generateBezierCoeff for quad, prism, hexa
-
-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.numFaces = -1;
-
- switch (tag){
- case MSH_PNT:
- //B.numFaces = 1;
- B.exposants = generate1DExposants(0);
- B.points = generate1DPoints(0);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 0);
- break;
- case MSH_LIN_2 :
- //B.numFaces = 2;
- B.exposants = generate1DExposants(0);
- B.points = generate1DPoints(0);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 0);
- //generate1dVertexClosure(B.closures);
- break;
- case MSH_LIN_3 :
- //B.numFaces = 2;
- B.exposants = generate1DExposants(1);
- B.points = generate1DPoints(1);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 1);
- //generate1dVertexClosure(B.closures);
- break;
- case MSH_LIN_4:
- //B.numFaces = 2;
- B.exposants = generate1DExposants(2);
- B.points = generate1DPoints(2);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 2);
- //generate1dVertexClosure(B.closures);
- break;
- case MSH_LIN_5:
- //B.numFaces = 2;
- B.exposants = generate1DExposants(3);
- B.points = generate1DPoints(3);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 3);
- //generate1dVertexClosure(B.closures);
- break;
- case MSH_LIN_6:
- //B.numFaces = 2;
- B.exposants = generate1DExposants(4);
- B.points = generate1DPoints(4);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 4);
- //generate1dVertexClosure(B.closures);
- break;
- case MSH_LIN_7:
- //B.numFaces = 2;
- B.exposants = generate1DExposants(5);
- B.points = generate1DPoints(5);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 5);
- //generate1dVertexClosure(B.closures);
- break;
- case MSH_LIN_8:
- //B.numFaces = 2;
- B.exposants = generate1DExposants(6);
- B.points = generate1DPoints(6);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 6);
- //generate1dVertexClosure(B.closures);
- break;
- case MSH_LIN_9:
- //B.numFaces = 2;
- B.exposants = generate1DExposants(7);
- B.points = generate1DPoints(7);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 7);
- //generate1dVertexClosure(B.closures);
- break;
- case MSH_LIN_10:
- //B.numFaces = 2;
- B.exposants = generate1DExposants(8);
- B.points = generate1DPoints(8);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 8);
- //generate1dVertexClosure(B.closures);
- break;
- case MSH_LIN_11:
- //B.numFaces = 2;
- B.exposants = generate1DExposants(9);
- B.points = generate1DPoints(9);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 9);
- //generate1dVertexClosure(B.closures);
- break;
- case MSH_TRI_3 :
- //B.numFaces = 3;
- B.exposants = generateExposantsTriangle(0);
- B.points = gmshGeneratePointsTriangle(0, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 0);
- //generate2dEdgeClosure(B.closures, 1);
- break;
- case MSH_TRI_6 :
- //B.numFaces = 3;
- B.exposants = generateExposantsTriangle(2);
- B.points = gmshGeneratePointsTriangle(2, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 2);
- //generate2dEdgeClosure(B.closures, 2);
- break;
- case MSH_TRI_9 :
- //B.numFaces = 3;
- B.exposants = generateExposantsTriangle(4);
- B.points = gmshGeneratePointsTriangle(4, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 4);
- //generate2dEdgeClosure(B.closures, 3);
- break;
- case MSH_TRI_10 :
- //B.numFaces = 3;
- B.exposants = generateExposantsTriangle(4);
- B.points = gmshGeneratePointsTriangle(4, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 4);
- //generate2dEdgeClosure(B.closures, 3);
- break;
- case MSH_TRI_12 :
- //B.numFaces = 3;
- B.exposants = generateExposantsTriangle(6);
- B.points = gmshGeneratePointsTriangle(6, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 6);
- //generate2dEdgeClosure(B.closures, 4);
- break;
- case MSH_TRI_15 :
- //B.numFaces = 3;
- B.exposants = generateExposantsTriangle(6);
- B.points = gmshGeneratePointsTriangle(6, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 6);
- //generate2dEdgeClosure(B.closures, 4);
- break;
- case MSH_TRI_15I :
- //B.numFaces = 3;
- B.exposants = generateExposantsTriangle(8);
- B.points = gmshGeneratePointsTriangle(8, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 8);
- //generate2dEdgeClosure(B.closures, 5);
- break;
- case MSH_TRI_21 :
- //B.numFaces = 3;
- B.exposants = generateExposantsTriangle(8);
- B.points = gmshGeneratePointsTriangle(8, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 8);
- //generate2dEdgeClosure(B.closures, 5);
- break;
- case MSH_TRI_28 :
- //B.numFaces = 3;
- B.exposants = generateExposantsTriangle(10);
- B.points = gmshGeneratePointsTriangle(10, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 10);
- //generate2dEdgeClosure(B.closures, 6);
- break;
- case MSH_TRI_36 :
- //B.numFaces = 3;
- B.exposants = generateExposantsTriangle(12);
- B.points = gmshGeneratePointsTriangle(12, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 12);
- //generate2dEdgeClosure(B.closures, 7);
- break;
- case MSH_TRI_45 :
- //B.numFaces = 3;
- B.exposants = generateExposantsTriangle(14);
- B.points = gmshGeneratePointsTriangle(14, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 14);
- //generate2dEdgeClosure(B.closures, 8);
- break;
- case MSH_TRI_55 :
- //B.numFaces = 3;
- B.exposants = generateExposantsTriangle(16);
- B.points = gmshGeneratePointsTriangle(16, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 16);
- //generate2dEdgeClosure(B.closures, 9);
- break;
- case MSH_TRI_66 :
- //B.numFaces = 3;
- B.exposants = generateExposantsTriangle(18);
- B.points = gmshGeneratePointsTriangle(18, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 18);
- //generate2dEdgeClosure(B.closures, 10);
- break;
- case MSH_TRI_18 :
- //B.numFaces = 3;
- B.exposants = generateExposantsTriangle(10);
- B.points = gmshGeneratePointsTriangle(10, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 10);
- //generate2dEdgeClosure(B.closures, 6);
- break;
- case MSH_TRI_21I :
- //B.numFaces = 3;
- B.exposants = generateExposantsTriangle(12);
- B.points = gmshGeneratePointsTriangle(12, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 12);
- //generate2dEdgeClosure(B.closures, 7);
- break;
- case MSH_TRI_24 :
- //B.numFaces = 3;
- B.exposants = generateExposantsTriangle(14);
- B.points = gmshGeneratePointsTriangle(14, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 14);
- //generate2dEdgeClosure(B.closures, 8);
- break;
- case MSH_TRI_27 :
- //B.numFaces = 3;
- B.exposants = generateExposantsTriangle(16);
- B.points = gmshGeneratePointsTriangle(16, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 16);
- //generate2dEdgeClosure(B.closures, 9);
- break;
- case MSH_TRI_30 :
- //B.numFaces = 3;
- B.exposants = generateExposantsTriangle(18);
- B.points = gmshGeneratePointsTriangle(18, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 18);
- //generate2dEdgeClosure(B.closures, 10);
- break;
- case MSH_TET_4 :
- //B.numFaces = 4;
- B.exposants = generateExposantsTetrahedron(0);
- B.points = gmshGeneratePointsTetrahedron(0, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 0);
- //generate3dFaceClosure(B.closures, 1);
- break;
- case MSH_TET_10 :
- //B.numFaces = 4;
- B.exposants = generateExposantsTetrahedron(3);
- B.points = gmshGeneratePointsTetrahedron(3, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 3);
- //generate3dFaceClosure(B.closures, 2);
- break;
- case MSH_TET_20 :
- //B.numFaces = 4;
- B.exposants = generateExposantsTetrahedron(6);
- B.points = gmshGeneratePointsTetrahedron(6, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 6);
- //generate3dFaceClosure(B.closures, 3);
- break;
- case MSH_TET_35 :
- //B.numFaces = 4;
- B.exposants = generateExposantsTetrahedron(9);
- B.points = gmshGeneratePointsTetrahedron(9, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 9);
- //generate3dFaceClosure(B.closures, 4);
- break;
- case MSH_TET_34 :
- //B.numFaces = 4;
- B.exposants = generateExposantsTetrahedron(9);
- B.points = gmshGeneratePointsTetrahedron(9, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 9);
- //generate3dFaceClosure(B.closures, 4);
- break;
- case MSH_TET_52 :
- //B.numFaces = 4;
- B.exposants = generateExposantsTetrahedron(12);
- B.points = gmshGeneratePointsTetrahedron(12, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 12);
- //generate3dFaceClosure(B.closures, 5);
- break;
- case MSH_TET_56 :
- //B.numFaces = 4;
- B.exposants = generateExposantsTetrahedron(12);
- B.points = gmshGeneratePointsTetrahedron(12, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 12);
- //generate3dFaceClosure(B.closures, 5);
- break;
- case MSH_TET_84 :
- //B.numFaces = 4;
- B.exposants = generateExposantsTetrahedron(15);
- B.points = gmshGeneratePointsTetrahedron(15, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 15);
- //generate3dFaceClosure(B.closures, 6);
- break;
- case MSH_TET_120 :
- //B.numFaces = 4;
- B.exposants = generateExposantsTetrahedron(18);
- B.points = gmshGeneratePointsTetrahedron(18, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 18);
- //generate3dFaceClosure(B.closures, 7);
- break;
- case MSH_TET_165 :
- //B.numFaces = 4;
- B.exposants = generateExposantsTetrahedron(21);
- B.points = gmshGeneratePointsTetrahedron(21, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 21);
- //generate3dFaceClosure(B.closures, 8);
- break;
- case MSH_TET_220 :
- //B.numFaces = 4;
- B.exposants = generateExposantsTetrahedron(24);
- B.points = gmshGeneratePointsTetrahedron(24, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 24);
- //generate3dFaceClosure(B.closures, 9);
- break;
- case MSH_TET_286 :
- //B.numFaces = 4;
- B.exposants = generateExposantsTetrahedron(27);
- B.points = gmshGeneratePointsTetrahedron(27, false);
- B.matrixLag2Bez = generateBezierFFCoefficientsSimplex(B.exposants, B.points, 27);
- //generate3dFaceClosure(B.closures, 10);
- break;
- case MSH_QUA_4 :
- //B.numFaces = 4;
- B.exposants = generateExposantsQuad(1);
- B.points = gmshGeneratePointsQuad(1,false);
- //generate2dEdgeClosure(B.closures, 1, 4);
- break;
- case MSH_QUA_9 :
- //B.numFaces = 4;
- B.exposants = generateExposantsQuad(3);
- B.points = gmshGeneratePointsQuad(3,false);
- //generate2dEdgeClosure(B.closures, 2, 4);
- break;
- case MSH_QUA_16 :
- //B.numFaces = 4;
- B.exposants = generateExposantsQuad(5);
- B.points = gmshGeneratePointsQuad(5,false);
- //generate2dEdgeClosure(B.closures, 3, 4);
- break;
- case MSH_QUA_25 :
- //B.numFaces = 4;
- B.exposants = generateExposantsQuad(7);
- B.points = gmshGeneratePointsQuad(7,false);
- //generate2dEdgeClosure(B.closures, 4, 4);
- break;
- case MSH_QUA_36 :
- //B.numFaces = 4;
- B.exposants = generateExposantsQuad(9);
- B.points = gmshGeneratePointsQuad(9,false);
- //generate2dEdgeClosure(B.closures, 5, 4);
- break;
- case MSH_QUA_49 :
- //B.numFaces = 4;
- B.exposants = generateExposantsQuad(11);
- B.points = gmshGeneratePointsQuad(11,false);
- //generate2dEdgeClosure(B.closures, 6, 4);
- break;
- case MSH_QUA_64 :
- //B.numFaces = 4;
- B.exposants = generateExposantsQuad(13);
- B.points = gmshGeneratePointsQuad(13,false);
- //generate2dEdgeClosure(B.closures, 7, 4);
- break;
- case MSH_QUA_81 :
- //B.numFaces = 4;
- B.exposants = generateExposantsQuad(15);
- B.points = gmshGeneratePointsQuad(15,false);
- //generate2dEdgeClosure(B.closures, 8, 4);
- break;
- case MSH_QUA_100 :
- //B.numFaces = 4;
- B.exposants = generateExposantsQuad(17);
- B.points = gmshGeneratePointsQuad(17,false);
- //generate2dEdgeClosure(B.closures, 9, 4);
- break;
- case MSH_QUA_121 :
- //B.numFaces = 4;
- B.exposants = generateExposantsQuad(19);
- B.points = gmshGeneratePointsQuad(19,false);
- //generate2dEdgeClosure(B.closures, 10, 4);
- break;
- case MSH_QUA_8 :
- //B.numFaces = 4;
- B.exposants = generateExposantsQuad(3);
- B.points = gmshGeneratePointsQuad(3,false);
- //generate2dEdgeClosure(B.closures, 2, 4);
- break;
- case MSH_QUA_12 :
- //B.numFaces = 4;
- B.exposants = generateExposantsQuad(5);
- B.points = gmshGeneratePointsQuad(5,false);
- //generate2dEdgeClosure(B.closures, 3, 4);
- break;
- case MSH_QUA_16I :
- //B.numFaces = 4;
- B.exposants = generateExposantsQuad(7);
- B.points = gmshGeneratePointsQuad(7,false);
- //generate2dEdgeClosure(B.closures, 4, 4);
- break;
- case MSH_QUA_20 :
- //B.numFaces = 4;
- B.exposants = generateExposantsQuad(9);
- B.points = gmshGeneratePointsQuad(9,false);
- //generate2dEdgeClosure(B.closures, 5, 4);
- break;
- case MSH_QUA_24 :
- //B.numFaces = 4;
- B.exposants = generateExposantsQuad(11);
- B.points = gmshGeneratePointsQuad(11,false);
- //generate2dEdgeClosure(B.closures, 6, 4);
- break;
- case MSH_QUA_28 :
- //B.numFaces = 4;
- B.exposants = generateExposantsQuad(13);
- B.points = gmshGeneratePointsQuad(13,false);
- //generate2dEdgeClosure(B.closures, 7, 4);
- break;
- case MSH_QUA_32 :
- //B.numFaces = 4;
- B.exposants = generateExposantsQuad(15);
- B.points = gmshGeneratePointsQuad(15,false);
- //generate2dEdgeClosure(B.closures, 8, 4);
- break;
- case MSH_QUA_36I :
- //B.numFaces = 4;
- B.exposants = generateExposantsQuad(17);
- B.points = gmshGeneratePointsQuad(17,false);
- //generate2dEdgeClosure(B.closures, 9, 4);
- break;
- case MSH_QUA_40 :
- //B.numFaces = 4;
- B.exposants = generateExposantsQuad(19);
- B.points = gmshGeneratePointsQuad(19,false);
- //generate2dEdgeClosure(B.closures, 10, 4);
- break;
- case MSH_PRI_6 : // first order
- //B.numFaces = 5;
- B.exposants = generateExposantsPrism(1);
- B.points = gmshGeneratePointsPrism(1, false);
- //generate3dFaceClosurePrism(B.closures, 1);
- break;
- case MSH_PRI_18 : // second order
- //B.numFaces = 5;
- B.exposants = generateExposantsPrism(4);
- B.points = gmshGeneratePointsPrism(4, false);
- //generate3dFaceClosurePrism(B.closures, 2);
- break;
-
- default :
- Msg::Error("Unknown function space %d: reverting to TET_4", tag);
- //B.numFaces = 4;
- B.exposants = generateExposantsTetrahedron(0);
- B.points = gmshGeneratePointsTetrahedron(0, false);
- //generate3dFaceClosure(B.closures, 1);
- break;
- }
-
- //B.matrixLag2Bez = generateLagrangeMonomialCoefficients(B.monomials, B.points);
-
- fs.insert(std::make_pair(tag, B));
- return &fs[tag];
-}
-//
-//
-//std::map<std::pair<int, int>, fullMatrix<double> > polynomialBases::injector;
-//
-//const fullMatrix<double> &polynomialBases::findInjector(int tag1, int tag2)
-//{
-// std::pair<int,int> key(tag1,tag2);
-// std::map<std::pair<int, int>, fullMatrix<double> >::const_iterator it = injector.find(key);
-// if (it != injector.end()) return it->second;
-//
-// const polynomialBasis& fs1 = *find(tag1);
-// const polynomialBasis& fs2 = *find(tag2);
-//
-// fullMatrix<double> inj(fs1.points.size1(), fs2.points.size1());
-//
-// double sf[256];
-//
-// for (int i = 0; i < fs1.points.size1(); i++) {
-// fs2.f(fs1.points(i, 0), fs1.points(i, 1), fs1.points(i, 2), sf);
-// for (int j = 0; j < fs2.points.size1(); j++) inj(i, j) = sf[j];
-// }
-//
-// injector.insert(std::make_pair(key, inj));
-// return injector[key];
-//}
-//
-//#include "Bindings.h"
-//void polynomialBasis::registerBindings(binding *b) {
-// classBinding *cb = b->addClass<polynomialBasis>("polynomialBasis");
-// cb->setDescription("polynomial shape functions for elements");
-// methodBinding *mb = cb->addMethod("f",(void (polynomialBasis::*)(fullMatrix<double>&, fullMatrix<double>&))&polynomialBasis::f);
-// mb->setDescription("evaluate the shape functions");
-// mb->setArgNames("nodes","values",NULL);
-// mb = cb->addMethod("find",&polynomialBases::find);
-// mb->setDescription("return the polynomial basis corresponding to an element type");
-// mb->setArgNames("elementType",NULL);
-//}
\ No newline at end of file
+// Gmsh - Copyright (C) 1997-2010 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];
+}
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index 885cf3b1a1e3e4e6fdb3dee158af83b6cb2d6b79..97924d136c6b0e1008a9e6e56309a3c5789ee60c 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -6,21 +6,22 @@
#ifndef _JACOBIAN_BASIS_H_
#define _JACOBIAN_BASIS_H_
-#include <math.h>
#include <map>
#include <vector>
#include "fullMatrix.h"
-#include <iostream>
class JacobianBasis
{
- public:
- fullMatrix<double> exposants; //exposants of Bezier FF
- fullMatrix<double> points; //sampling point
- fullMatrix<double> matrixLag2Bez;
- fullMatrix<double> gradShapeMatX;
- fullMatrix<double> gradShapeMatY;
- fullMatrix<double> gradShapeMatZ;
+ 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
@@ -31,207 +32,4 @@ class JacobianBases
static const JacobianBasis *find(int);
};
-//// presently those function spaces are only for simplices and quads;
-//// should be extended to other elements like hexes
-//class polynomialBasis
-//{
-// public:
-// typedef std::vector<std::vector<int> > clCont;
-// clCont closures;
-// fullMatrix<double> points;
-// fullMatrix<double> monomials;
-// fullMatrix<double> coefficients;
-// int numFaces;
-// // for a given face/edge, with both a sign and a rotation,
-// // give an ordered list of nodes on this face/edge
-// inline const std::vector<int> &getClosure(int id) const // return the closure of dimension dim
-// {
-// return closures[id];
-// }
-// inline int getClosureId(int iEl, int iSign=1, int iRot=0) const {
-// return iEl + numFaces*(iSign == 1 ? 0 : 1) + 2*numFaces*iRot;
-// }
-// inline void evaluateMonomials(double u, double v, double w, double p[]) const
-// {
-// for (int j = 0; j < monomials.size1(); j++) {
-// p[j] = pow_int(u, (int)monomials(j, 0));
-// if (monomials.size2() > 1) p[j] *= pow_int(v, (int)monomials(j, 1));
-// if (monomials.size2() > 2) p[j] *= pow_int(w, (int)monomials(j, 2));
-// }
-// }
-// inline void f(double u, double v, double w, double *sf) const
-// {
-// double p[256];
-// evaluateMonomials(u, v, w, p);
-// for (int i = 0; i < coefficients.size1(); i++) {
-// sf[i] = 0;
-// for (int j = 0; j < coefficients.size2(); j++) {
-// sf[i] += coefficients(i, j) * p[j];
-// }
-// }
-// }
-// // I would favour this interface that allows optimizations (not implemented) and is easier to bind
-// inline void f(fullMatrix<double> &coord, fullMatrix<double> &sf) {
-// double p[256];
-// sf.resize (coord.size1(), coefficients.size1());
-// for (int iPoint=0; iPoint< coord.size1(); iPoint++) {
-// evaluateMonomials(coord(iPoint,0), coord(iPoint,1), coord(iPoint,2), p);
-// for (int i = 0; i < coefficients.size1(); i++)
-// for (int j = 0; j < coefficients.size2(); j++)
-// sf(iPoint,i) += coefficients(i, j) * p[j];
-// }
-// }
-// inline void df(double u, double v, double w, double grads[][3]) const
-// {
-// switch (monomials.size2()) {
-// case 1:
-// for (int i = 0; i < coefficients.size1(); i++){
-// grads[i][0] = 0;
-// grads[i][1] = 0;
-// grads[i][2] = 0;
-// for(int j = 0; j < coefficients.size2(); j++){
-// if (monomials(j, 0) > 0)
-// grads[i][0] += coefficients(i, j) *
-// pow_int(u, (int)(monomials(j, 0) - 1)) * monomials(j, 0);
-// }
-// }
-// break;
-// case 2:
-// for (int i = 0; i < coefficients.size1(); i++){
-// grads[i][0] = 0;
-// grads[i][1] = 0;
-// grads[i][2] = 0;
-// for(int j = 0; j < coefficients.size2(); j++){
-// if (monomials(j, 0) > 0)
-// grads[i][0] += coefficients(i, j) *
-// pow_int(u, (int)(monomials(j, 0) - 1)) * monomials(j, 0) *
-// pow_int(v, (int)monomials(j, 1));
-// if (monomials(j, 1) > 0)
-// grads[i][1] += coefficients(i, j) *
-// pow_int(u, (int)monomials(j, 0)) *
-// pow_int(v, (int)(monomials(j, 1) - 1)) * monomials(j, 1);
-// }
-// }
-// break;
-// case 3:
-// for (int i = 0; i < coefficients.size1(); i++){
-// grads[i][0] = 0;
-// grads[i][1] = 0;
-// grads[i][2] = 0;
-// for(int j = 0; j < coefficients.size2(); j++){
-// if (monomials(j, 0) > 0)
-// grads[i][0] += coefficients(i, j) *
-// pow_int(u, (int)(monomials(j, 0) - 1)) * monomials(j, 0) *
-// pow_int(v, (int)monomials(j, 1)) *
-// pow_int(w, (int)monomials(j, 2));
-// if (monomials(j, 1) > 0)
-// grads[i][1] += coefficients(i, j) *
-// pow_int(u, (int)monomials(j, 0)) *
-// pow_int(v, (int)(monomials(j, 1) - 1)) * monomials(j, 1) *
-// pow_int(w, (int)monomials(j, 2));
-// if (monomials(j, 2) > 0)
-// grads[i][2] += coefficients(i, j) *
-// pow_int(u, (int)monomials(j, 0)) *
-// pow_int(v, (int)monomials(j, 1)) *
-// pow_int(w, (int)(monomials(j, 2) - 1)) * monomials(j, 2);
-// }
-// }
-// break;
-// }
-// }
-// inline void ddf(double u, double v, double w, double hess[][3][3]) const
-// {
-// switch (monomials.size2()) {
-// case 1:
-// for (int i = 0; i < coefficients.size1(); i++){
-// hess[i][0][0] = hess[i][0][1] = hess[i][0][2] = 0;
-// hess[i][1][0] = hess[i][1][1] = hess[i][1][2] = 0;
-// hess[i][2][0] = hess[i][2][1] = hess[i][2][2] = 0;
-//
-// for(int j = 0; j < coefficients.size2(); j++){
-// if (monomials(j, 0) > 1) // second derivative !=0
-// hess[i][0][0] += coefficients(i, j) * pow_int(u, (int)(monomials(j, 0) - 2)) *
-// monomials(j, 0) * (monomials(j, 0) - 1);
-// }
-// }
-// break;
-// case 2:
-// for (int i = 0; i < coefficients.size1(); i++){
-// hess[i][0][0] = hess[i][0][1] = hess[i][0][2] = 0;
-// hess[i][1][0] = hess[i][1][1] = hess[i][1][2] = 0;
-// hess[i][2][0] = hess[i][2][1] = hess[i][2][2] = 0;
-// for(int j = 0; j < coefficients.size2(); j++){
-// if (monomials(j, 0) > 1) // second derivative !=0
-// hess[i][0][0] += coefficients(i, j) * pow_int(u, (int)(monomials(j, 0) - 2)) *
-// monomials(j, 0) * (monomials(j, 0) - 1) * pow_int(v, (int)monomials(j, 1));
-// if ((monomials(j, 1) > 0) && (monomials(j, 0) > 0))
-// hess[i][0][1] += coefficients(i, j) *
-// pow_int(u, (int)monomials(j, 0) - 1) * monomials(j, 0) *
-// pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1);
-// if (monomials(j, 1) > 1)
-// hess[i][1][1] += coefficients(i, j) *
-// pow_int(u, (int)monomials(j, 0)) *
-// pow_int(v, (int)monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1) - 1);
-// }
-// hess[i][1][0] = hess[i][0][1];
-// }
-// break;
-// case 3:
-// for (int i = 0; i < coefficients.size1(); i++){
-// hess[i][0][0] = hess[i][0][1] = hess[i][0][2] = 0;
-// hess[i][1][0] = hess[i][1][1] = hess[i][1][2] = 0;
-// hess[i][2][0] = hess[i][2][1] = hess[i][2][2] = 0;
-// for(int j = 0; j < coefficients.size2(); j++){
-// if (monomials(j, 0) > 1)
-// hess[i][0][0] += coefficients(i, j) *
-// pow_int(u, (int)monomials(j, 0) - 2) * monomials(j, 0) * (monomials(j, 0)-1) *
-// pow_int(v, (int)monomials(j, 1)) *
-// pow_int(w, (int)monomials(j, 2));
-//
-// if ((monomials(j, 0) > 0) && (monomials(j, 1) > 0))
-// hess[i][0][1] += coefficients(i, j) *
-// pow_int(u, (int)monomials(j, 0) - 1) * monomials(j, 0) *
-// pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1) *
-// pow_int(w, (int)monomials(j, 2));
-// if ((monomials(j, 0) > 0) && (monomials(j, 2) > 0))
-// hess[i][0][2] += coefficients(i, j) *
-// pow_int(u, (int)monomials(j, 0) - 1) * monomials(j, 0) *
-// pow_int(v, (int)monomials(j, 1)) *
-// pow_int(w, (int)monomials(j, 2) - 1) * monomials(j, 2);
-// if (monomials(j, 1) > 1)
-// hess[i][1][1] += coefficients(i, j) *
-// pow_int(u, (int)monomials(j, 0)) *
-// pow_int(v, (int)monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1)-1) *
-// pow_int(w, (int)monomials(j, 2));
-// if ((monomials(j, 1) > 0) && (monomials(j, 2) > 0))
-// hess[i][1][2] += coefficients(i, j) *
-// pow_int(u, (int)monomials(j, 0)) *
-// pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1) *
-// pow_int(w, (int)monomials(j, 2) - 1) * monomials(j, 2);
-// if (monomials(j, 2) > 1)
-// hess[i][2][2] += coefficients(i, j) *
-// pow_int(u, (int)monomials(j, 0)) *
-// pow_int(v, (int)monomials(j, 1)) *
-// pow_int(w, (int)monomials(j, 2) - 2) * monomials(j, 2) * (monomials(j, 2) - 1);
-// }
-// hess[i][1][0] = hess[i][0][1];
-// hess[i][2][0] = hess[i][0][2];
-// hess[i][2][1] = hess[i][1][2];
-// }
-// break;
-// }
-// }
-// static void registerBindings(binding *b);
-//};
-//
-//class polynomialBases
-//{
-// private:
-// static std::map<int, polynomialBasis> fs;
-// static std::map<std::pair<int, int>, fullMatrix<double> > injector;
-// public :
-// static const polynomialBasis *find(int);
-// static const fullMatrix<double> &findInjector(int, int);
-//};
-
#endif
diff --git a/Numeric/fullMatrix.cpp b/Numeric/fullMatrix.cpp
index faa67bd3ec2aced3e2ca4369588598bd24f5e0b4..bfc871fa1e49ad627f71eea9454c5cd0a0fad88f 100644
--- a/Numeric/fullMatrix.cpp
+++ b/Numeric/fullMatrix.cpp
@@ -42,7 +42,7 @@ extern "C" {
}
template<>
-void fullVector<double>::axpy(fullVector<double> &x,double alpha)
+void fullVector<double>::axpy(const fullVector<double> &x,double alpha)
{
int M = _r, INCX = 1, INCY = 1;
F77NAME(daxpy)(&M, &alpha, x._data,&INCX, _data, &INCY);
@@ -220,7 +220,7 @@ bool fullMatrix<double>::luSolve(const fullVector<double> &rhs, fullVector<doubl
}
template<>
-bool fullMatrix<double>::invert(fullMatrix<double> &result)
+bool fullMatrix<double>::invert(fullMatrix<double> &result) const
{
int M = size1(), N = size2(), lda = size1(), info;
int *ipiv = new int[std::min(M, N)];
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index 184ed22064b6e82b84a864f87c8b4e8f2c79a552..fb812896a04356ab3f9dd148bf655fb6e7086ec6 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -47,6 +47,13 @@ class fullVector
}
return false;
}
+ void setAsProxy(const fullVector<scalar> &original, int r_start, int r)
+ {
+ if(_data)
+ delete [] _data;
+ _r = r;
+ _data = original._data + r_start;
+ }
inline scalar operator () (int i) const
{
return _data[i];
@@ -77,13 +84,17 @@ class fullVector
for(int i = 0; i < _r; ++i) s += _data[i] * other._data[i];
return s;
}
- void axpy(fullVector<scalar> &x, scalar alpha=1.)
+ void axpy(const fullVector<scalar> &x, scalar alpha=1.)
#if !defined(HAVE_BLAS)
{
for (int i = 0; i < _r; i++) _data[i] += alpha * x._data[i];
}
#endif
;
+ void multTByT(const fullVector<scalar> &x)
+ {
+ for (int i = 0; i < _r; i++) _data[i] *= x._data[i];
+ }
void print(const char *name="") const
{
printf("Printing vector %s:\n", name);
@@ -260,6 +271,16 @@ class fullMatrix
for(int j = j0, destj = destj0; j < j0 + nj; j++, destj++)
(*this)(desti, destj) = a(i, j);
}
+ void copy(const fullMatrix<scalar> &a)
+ {
+ if(_data && _own_data)
+ delete [] _data;
+ _r = a._r;
+ _c = a._c;
+ _data = new scalar[_r * _c];
+ _own_data = true;
+ setAll(a);
+ }
void mult_naive(const fullMatrix<scalar> &b, fullMatrix<scalar> &c)const
{
c.scale(0.);
@@ -276,6 +297,10 @@ class fullMatrix
}
#endif
;
+ void multTByT(const fullMatrix<scalar> &a)
+ {
+ for (int i = 0; i < _r * _c; i++) _data[i] *= a._data[i];
+ }
void gemm_naive(const fullMatrix<scalar> &a, const fullMatrix<scalar> &b,
scalar alpha=1., scalar beta=1.)
{
@@ -386,7 +411,7 @@ class fullMatrix
}
#endif
;
- bool invert(fullMatrix<scalar> &result)
+ bool invert(fullMatrix<scalar> &result) const
#if !defined(HAVE_LAPACK)
{
Msg::Error("LU factorization requires LAPACK");
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 474f4f025098e6db7fc1d0dbe20fbddeb56f2d51..6fc0e592356a92ff9a375986524f08c6d8a1c36c 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -730,8 +730,8 @@ static fullMatrix<double> generateLagrangeMonomialCoefficients
{
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() );
+ monomial.size1(),point.size1(),
+ monomial.size2(),point.size2() );
return fullMatrix<double>(1, 1);
}
diff --git a/Plugin/AnalyseCurvedMesh.cpp b/Plugin/AnalyseCurvedMesh.cpp
index 0f0d3ba6ec6dc1f7b3600a0f15b227cf7d2b2c48..3edb3352987226db9f4d21d7cafec4590d4e4c6e 100644
--- a/Plugin/AnalyseCurvedMesh.cpp
+++ b/Plugin/AnalyseCurvedMesh.cpp
@@ -2,19 +2,21 @@
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
-//
-// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
-
-#include "GmshDefines.h"
-#include <stdlib.h>
+
#include "Gmsh.h"
-#include "MTriangle.h"
-#include "GmshConfig.h"
-#include "GModel.h"
-#include "polynomialBasis.h"
+#include "GModel.h"
+#include "GmshMessage.h"
#include "JacobianBasis.h"
+#include "polynomialBasis.h"
#include "AnalyseCurvedMesh.h"
-#include "fullMatrix.h"
+#define UNDEF_JAC_TAG -999
+
+
+StringXNumber JacobianOptions_Number[] = {
+ {GMSH_FULLRC, "Method", NULL, 1},
+ {GMSH_FULLRC, "MaxDepth", NULL, 5}
+};
+
extern "C"
{
@@ -22,12 +24,38 @@ extern "C"
{
return new GMSH_AnalyseCurvedMeshPlugin();
}
-}
+}
+int GMSH_AnalyseCurvedMeshPlugin::getNbOptions() const
+{
+ return sizeof(JacobianOptions_Number) / sizeof(StringXNumber);
+}
+
+StringXNumber *GMSH_AnalyseCurvedMeshPlugin::getOption(int iopt)
+{
+ return &JacobianOptions_Number[iopt];
+}
std::string GMSH_AnalyseCurvedMeshPlugin::getHelp() const
{
- return "Plugin(AnalyseCurvedMesh) check the jacobian of all elements of the greater dimension. Elements for which we are absolutely certain that they are good (positive jacobian) are ignored. Others are only suppose to be wrong."
- "Plugin(AnalyseCurvedMesh) write in the console which elements are supposed to be wrong. (if labels of analysed type of elements are set visible, only supposed wrong elements will be visible)";
+ return "Plugin(AnalyseCurvedMesh) check the jacobian of all elements of the greater dimension. "
+ "Elements for which we are absolutely certain that they are good (positive jacobian) are ignored. "
+ "Others are wrong or suppose to be wrong.\n\n"
+ "Plugin(AnalyseCurvedMesh) write in the console which elements are wrong. "
+ "(if labels of analysed type of elements are set visible, only wrong elements will be visible)\n\n"
+ "method = 1 or 2\n"
+ "maxDepth >= 0 (> 1 is better)\n\n"
+ "maxDepth = 0 : only sampling of the jacobian\n"
+ "maxDepth = 1 : also calculate control points\n"
+ "maxDepth = 2+ : also decompose control points\n";
+}
+
+// Miscellaneous method
+int max(const std::vector<int> &vec)
+{
+ int max = vec[0];
+ for (unsigned int i = 1; i < vec.size(); i++)
+ if (vec[i] > max) max = vec[i];
+ return max;
}
static double computeDeterminant(MElement *el, double jac[3][3])
@@ -43,8 +71,8 @@ static double computeDeterminant(MElement *el, double jac[3][3])
return jac[0][0] * jac[1][1] * jac[2][2] + jac[0][2] * jac[1][0] * jac[2][1] +
jac[0][1] * jac[1][2] * jac[2][0] - jac[0][2] * jac[1][1] * jac[2][0] -
jac[0][0] * jac[1][2] * jac[2][1] - jac[0][1] * jac[1][0] * jac[2][2];
- default:
- return 1;
+ default:
+ return 1;
}
}
@@ -67,219 +95,706 @@ double getJacobian(double gsf[][3], double jac[3][3], MElement *el)
return computeDeterminant(el, jac);
}
+void setJacobian(MElement *el, const JacobianBasis *jfs, fullVector<double> &jacobian)
+{
+ int numVertices = el->getNumVertices();
+ fullVector<double> nodesX(numVertices);
+ fullVector<double> nodesY;
+ fullVector<double> nodesZ;
+ fullVector<double> interm1;
+ fullVector<double> interm2;
+
+ switch (el->getDim()) {
+
+ case 1 :
+ for (int i = 0; i < numVertices; i++) {
+ nodesX(i) = el->getVertex(i)->x();
+ }
+ jfs->gradShapeMatX.mult(nodesX, jacobian);
+ break;
+
+
+ case 2 :
+
+ nodesY.resize(numVertices);
+ interm1.resize(jacobian.size());
+ interm2.resize(jacobian.size());
+
+ for (int i = 0; i < numVertices; i++) {
+ nodesX(i) = el->getVertex(i)->x();
+ nodesY(i) = el->getVertex(i)->y();
+ }
+
+ jfs->gradShapeMatX.mult(nodesX, jacobian);
+ jfs->gradShapeMatY.mult(nodesY, interm2);
+ jacobian.multTByT(interm2);
+
+ jfs->gradShapeMatY.mult(nodesX, interm1);
+ jfs->gradShapeMatX.mult(nodesY, interm2);
+ interm1.multTByT(interm2);
+
+ jacobian.axpy(interm1, -1);
+ break;
+
+
+ case 3 :
+
+ nodesY.resize(numVertices);
+ nodesZ.resize(numVertices);
+ interm1.resize(jacobian.size());
+ interm2.resize(jacobian.size());
+
+ for (int i = 0; i < numVertices; i++) {
+ nodesX(i) = el->getVertex(i)->x();
+ nodesY(i) = el->getVertex(i)->y();
+ nodesZ(i) = el->getVertex(i)->z();
+ }
+
+ jfs->gradShapeMatX.mult(nodesX, jacobian);
+ jfs->gradShapeMatY.mult(nodesY, interm2);
+ jacobian.multTByT(interm2);
+ jfs->gradShapeMatZ.mult(nodesZ, interm2);
+ jacobian.multTByT(interm2);
+
+ jfs->gradShapeMatX.mult(nodesY, interm1);
+ jfs->gradShapeMatY.mult(nodesZ, interm2);
+ interm1.multTByT(interm2);
+ jfs->gradShapeMatZ.mult(nodesX, interm2);
+ interm1.multTByT(interm2);
+ jacobian.axpy(interm1, 1);
+
+ jfs->gradShapeMatX.mult(nodesZ, interm1);
+ jfs->gradShapeMatY.mult(nodesX, interm2);
+ interm1.multTByT(interm2);
+ jfs->gradShapeMatZ.mult(nodesY, interm2);
+ interm1.multTByT(interm2);
+ jacobian.axpy(interm1, 1);
+
+
+ jfs->gradShapeMatX.mult(nodesY, interm1);
+ jfs->gradShapeMatY.mult(nodesX, interm2);
+ interm1.multTByT(interm2);
+ jfs->gradShapeMatZ.mult(nodesZ, interm2);
+ interm1.multTByT(interm2);
+ jacobian.axpy(interm1, -1);
+
+ jfs->gradShapeMatX.mult(nodesZ, interm1);
+ jfs->gradShapeMatY.mult(nodesY, interm2);
+ interm1.multTByT(interm2);
+ jfs->gradShapeMatZ.mult(nodesX, interm2);
+ interm1.multTByT(interm2);
+ jacobian.axpy(interm1, -1);
+
+ jfs->gradShapeMatX.mult(nodesX, interm1);
+ jfs->gradShapeMatY.mult(nodesZ, interm2);
+ interm1.multTByT(interm2);
+ jfs->gradShapeMatZ.mult(nodesY, interm2);
+ interm1.multTByT(interm2);
+ jacobian.axpy(interm1, -1);
+ }
+}
+
+
+void setJacobian(MElement *const *el, const JacobianBasis *jfs, fullMatrix<double> &jacobian)
+{
+ int numEl = jacobian.size2();
+ int numVertices = el[0]->getNumVertices();
+ fullMatrix<double> nodesX(numVertices,numEl);
+ fullMatrix<double> nodesY;
+ fullMatrix<double> nodesZ;
+ fullMatrix<double> interm1;
+ fullMatrix<double> interm2;
+
+ switch (el[0]->getDim()) {
+
+ case 1 :
+ for (int j = 0; j < numEl; j++) {
+ for (int i = 0; i < numVertices; i++) {
+ nodesX(i,j) = el[j]->getVertex(i)->x();
+ }
+ }
+ jfs->gradShapeMatX.mult(nodesX, jacobian);
+ break;
+
+
+ case 2 :
+
+ nodesY.resize(numVertices,numEl);
+ interm1.resize(jacobian.size1(),jacobian.size2());
+ interm2.resize(jacobian.size1(),jacobian.size2());
+
+ for (int j = 0; j < numEl; j++) {
+ for (int i = 0; i < numVertices; i++) {
+ nodesX(i,j) = el[j]->getVertex(i)->x();
+ nodesY(i,j) = el[j]->getVertex(i)->y();
+ }
+ }
+
+ jfs->gradShapeMatX.mult(nodesX, jacobian);
+ jfs->gradShapeMatY.mult(nodesY, interm2);
+ jacobian.multTByT(interm2);
+
+ jfs->gradShapeMatY.mult(nodesX, interm1);
+ jfs->gradShapeMatX.mult(nodesY, interm2);
+ interm1.multTByT(interm2);
+
+ jacobian.add(interm1, -1);
+ break;
+
+
+ case 3 :
+
+ nodesY.resize(numVertices,numEl);
+ nodesZ.resize(numVertices,numEl);
+ interm1.resize(jacobian.size1(),jacobian.size2());
+ interm2.resize(jacobian.size1(),jacobian.size2());
+
+ for (int j = 0; j < numEl; j++) {
+ for (int i = 0; i < numVertices; i++) {
+ nodesX(i,j) = el[j]->getVertex(i)->x();
+ nodesY(i,j) = el[j]->getVertex(i)->y();
+ nodesZ(i,j) = el[j]->getVertex(i)->z();
+ }
+ }
+
+ jfs->gradShapeMatX.mult(nodesX, jacobian);
+ jfs->gradShapeMatY.mult(nodesY, interm2);
+ jacobian.multTByT(interm2);
+ jfs->gradShapeMatZ.mult(nodesZ, interm2);
+ jacobian.multTByT(interm2);
+
+ jfs->gradShapeMatX.mult(nodesY, interm1);
+ jfs->gradShapeMatY.mult(nodesZ, interm2);
+ interm1.multTByT(interm2);
+ jfs->gradShapeMatZ.mult(nodesX, interm2);
+ interm1.multTByT(interm2);
+ jacobian.add(interm1, 1);
+
+ jfs->gradShapeMatX.mult(nodesZ, interm1);
+ jfs->gradShapeMatY.mult(nodesX, interm2);
+ interm1.multTByT(interm2);
+ jfs->gradShapeMatZ.mult(nodesY, interm2);
+ interm1.multTByT(interm2);
+ jacobian.add(interm1, 1);
+
+
+ jfs->gradShapeMatX.mult(nodesY, interm1);
+ jfs->gradShapeMatY.mult(nodesX, interm2);
+ interm1.multTByT(interm2);
+ jfs->gradShapeMatZ.mult(nodesZ, interm2);
+ interm1.multTByT(interm2);
+ jacobian.add(interm1, -1);
+
+ jfs->gradShapeMatX.mult(nodesZ, interm1);
+ jfs->gradShapeMatY.mult(nodesY, interm2);
+ interm1.multTByT(interm2);
+ jfs->gradShapeMatZ.mult(nodesX, interm2);
+ interm1.multTByT(interm2);
+ jacobian.add(interm1, -1);
+
+ jfs->gradShapeMatX.mult(nodesX, interm1);
+ jfs->gradShapeMatY.mult(nodesZ, interm2);
+ interm1.multTByT(interm2);
+ jfs->gradShapeMatZ.mult(nodesY, interm2);
+ interm1.multTByT(interm2);
+ jacobian.add(interm1, -1);
+ }
+}
+
+
+// Execution
PView *GMSH_AnalyseCurvedMeshPlugin::execute(PView *v)
{
- Msg::Info("AnalyseCurvedMeshPlugin : Starting analyse.");
- int numBadEl = 0;
- int numAnalysedEl = 0;
-
- GModel *m = GModel::current();
-
- JacobianBases::find(MSH_QUA_8);
- JacobianBases::find(MSH_QUA_9);
- JacobianBases::find(MSH_TET_20);
- JacobianBases::find(MSH_PRI_18);
-
- switch (m->getDim()) {
-
- case 3 :
- {
- Msg::Info("Only 3D elements will be analyse.");
- for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); it++) {
- GRegion *r = *it;
-
- unsigned int numType[5] = {0, 0, 0, 0, 0};
- r->getNumMeshElements(numType);
-
- for(int type = 0; type < 5; type++) { // loop over all types of elements
- MElement *const *el = r->getStartElementType(type);
-
- for(unsigned int i = 0; i < numType[type]; i++) { // loop over all elements of a type
- numAnalysedEl++;
- if (checkJacobian((*(el+i)), 0) <= 0) numBadEl++;
- }
- }
- }
- break;
- }
-
- case 2 :
- {
- Msg::Info("Only 2D elements will be analyse.");
- Msg::Warning("2D elements must be in a z=cst plane ! If they aren't, results won't be correct.");
- for (GModel::fiter it = m->firstFace(); it != m->lastFace(); it++) {
- GFace *f = *it;
-
- unsigned int numType[3] = {0, 0, 0};
- f->getNumMeshElements(numType);
-
- for (int type = 0; type < 3; type++) {
- MElement *const *el = f->getStartElementType(type);
-
- for (unsigned int i = 0; i < numType[type]; i++) {
- numAnalysedEl++;
- if (checkJacobian((*(el+i)), 0) <= 0) numBadEl++;
- }
- }
- }
- break;
- }
-
- case 1 :
- {
- Msg::Info("Only 1D elements will be analyse.");
- Msg::Warning("1D elements must be on a y=cst & z=cst line ! If they aren't, results won't be correct.");
- for (GModel::eiter it = m->firstEdge(); it != m->lastEdge(); it++) {
- GEdge *e = *it;
-
- unsigned int *numElement = {0};
- e->getNumMeshElements(numElement);
-
- MElement *const *el = e->getStartElementType(0);
-
- for (unsigned int i = 0; i < *numElement; i++) {
- numAnalysedEl++;
- if (checkJacobian((*(el+i)), 0) <= 0) numBadEl++;
- }
- }
- break;
- }
-
- default :
- {
- Msg::Error("I can't analyse any element.");
- }
- }
-
- Msg::Info("%d elements have been analysed.", numAnalysedEl);
- Msg::Info("%d elements may be bad.", numBadEl);
- Msg::Info("AnalyseCurvedMeshPlugin : Job finished.");
+ Msg::Info("AnalyseCurvedMeshPlugin : Starting analyse.");
+ int numBadEl = 0;
+ int numUncertain = 0;
+ int numAnalysedEl = 0;
+ std::vector<int> tag;
+ int method = (int)JacobianOptions_Number[0].def;
+ int maxDepth = (int)JacobianOptions_Number[1].def;
+
+ if (method < 1 || method > 2) {
+#if defined(HAVE_BLAS)
+ method = 2;
+#else
+ method = 1;
+#endif
+ }
+
+ GModel *m = GModel::current();
+ switch (m->getDim()) {
+
+ case 3 :
+ Msg::Info("Only 3D elements will be analyse.");
+ for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); it++) {
+ GRegion *r = *it;
+
+ unsigned int numType[5] = {0, 0, 0, 0, 0};
+ r->getNumMeshElements(numType);
+
+ for(int type = 0; type < 5; type++) {
+ MElement *const *el = r->getStartElementType(type);
+
+ int *a;
+ a = checkJacobian(el, numType[type], maxDepth, method);
+ numUncertain += a[0];
+ numBadEl += a[1];
+ numAnalysedEl += numType[type];
+ delete[] a;
+
+ /*for(unsigned int i = 0; i < numType[type]; i++) {
+ numAnalysedEl++;
+ if (checkJacobian(el[i], maxDepth) <= 0) numBadEl++;
+ }*/
+ }
+ }
+ break;
+
+ case 2 :
+ Msg::Info("Only 2D elements will be analyse.");
+ Msg::Warning("2D elements must be in a z=cst plane ! If they aren't, results won't be correct.");
+ for (GModel::fiter it = m->firstFace(); it != m->lastFace(); it++) {
+ GFace *f = *it;
+
+ unsigned int numType[3] = {0, 0, 0};
+ f->getNumMeshElements(numType);
+
+ for (int type = 0; type < 3; type++) {
+ MElement *const *el = f->getStartElementType(type);
+
+ int *a;
+ a = checkJacobian(el, numType[type], maxDepth, method);
+ numUncertain += a[0];
+ numBadEl += a[1];
+ numAnalysedEl += numType[type];
+ delete[] a;
+
+ /*for (unsigned int i = 0; i < numType[type]; i++) {
+ numAnalysedEl++;
+ if (checkJacobian(el[i], maxDepth) <= 0) numBadEl++;
+ }*/
+ }
+ }
+ break;
+
+ case 1 :
+ Msg::Info("Only 1D elements will be analyse.");
+ Msg::Warning("1D elements must be on a y=cst & z=cst line ! If they aren't, results won't be correct.");
+ for (GModel::eiter it = m->firstEdge(); it != m->lastEdge(); it++) {
+ GEdge *e = *it;
+
+ unsigned int numElement = e->getNumMeshElements();
+ MElement *const *el = e->getStartElementType(0);
+
+ int *a;
+ a = checkJacobian(el, numElement, maxDepth, method);
+ numUncertain += a[0];
+ numBadEl += a[1];
+ numAnalysedEl += numElement;
+ delete[] a;
+
+ /*for (unsigned int i = 0; i < numElement; i++) {
+ numAnalysedEl++;
+ if (checkJacobian(el[i], maxDepth) <= 0) numBadEl++;
+ }*/
+ }
+ break;
+
+ default :
+ Msg::Error("I can't analyse any element.");
+ }
+
+
+ //Set all visibility of smaller dimension elements to 0
+ switch (m->getDim()) {
+ case 2 :
+ for (GModel::fiter it = m->firstFace(); it != m->lastFace(); it++) {
+ GFace *f = *it;
+
+ unsigned int numType[3] = {0, 0, 0};
+ f->getNumMeshElements(numType);
+
+ for (int type = 0; type < 3; type++) {
+ MElement *const *el = f->getStartElementType(type);
+ for (unsigned int i = 0; i < numType[type]; i++) {
+ el[i]->setVisibility(0);
+ }
+ }
+ }
+ case 1 :
+ for (GModel::eiter it = m->firstEdge(); it != m->lastEdge(); it++) {
+ GEdge *e = *it;
+
+ unsigned int numElement = e->getNumMeshElements();
+ MElement *const *el = e->getStartElementType(0);
+
+ for (unsigned int i = 0; i < numElement; i++) {
+ el[i]->setVisibility(0);
+ }
+ }
+ break;
+ }
+
+ Msg::Info("%d elements have been analysed.", numAnalysedEl);
+ Msg::Info("%d elements were bad.", numBadEl);
+ Msg::Info("%d elements were undetermined.", numUncertain);
+ Msg::Info("AnalyseCurvedMeshPlugin : Job finished.");
return 0;
}
-bool GMSH_AnalyseCurvedMeshPlugin::isJacPositive(MElement *el)
+int *GMSH_AnalyseCurvedMeshPlugin::checkJacobian
+(MElement *const *el, int numEl, int maxDepth, int method)
{
- const polynomialBasis *fs = el->getFunctionSpace(-1);
- if (!fs) {
- Msg::Error("Function space not implemented for type of element %d", el->getNum());
- return false;
- }
- const JacobianBasis *jfs = el->getJacobianFuncSpace(-1);
- if (!jfs) {
- Msg::Error("Jacobian function space not implemented for type of element %d", el->getNum());
- return false;
- }
- int numSamplingPt = jfs->points.size1();
- int dim = jfs->points.size2();
- bool isPositive = true;
- fullVector<double> jacobian(numSamplingPt);
-
- for (int i = 0; i < numSamplingPt; i++) {
- double gsf[256][3];
- switch (dim) {
- case 1:
- fs->df(jfs->points(i,0),0,0, gsf);
- break;
- case 2:
- fs->df(jfs->points(i,0),jfs->points(i,1),0, gsf);
- break;
- case 3:
- fs->df(jfs->points(i,0),jfs->points(i,1),jfs->points(i,2), gsf);
- break;
- default:
- Msg::Error("Can't get the gradient for %dD elements.", dim);
- return false;
- }
- double jac[3][3];
- jacobian(i) = getJacobian(gsf, jac, el);
- }
-
- fullVector<double> jacBez(jacobian.size());
-
- jfs->matrixLag2Bez.mult(jacobian, jacBez);
-
- for (int i = 0; i < jacBez.size(); i++) {
- if (jacBez(i) < 0.) isPositive = false;
- }
-
- //Msg::Info("%d", el->getNum());
- //if (el->getNum() == 582 || el->getNum() == 584) {
- // for (int i = 0; i < jfs->points.size1(); i++) {
- // Msg::Info("jacBez[%d] = %lf, jacobian[%d] = %lf", i, jacBez(i), i, jacobian(i));
- // }
- //}
-
- return isPositive;
+ int *a = new int[2];
+ a[0] = a[1] = 0;
+ if (numEl <= 0) return a;
+
+ switch (method) {
+
+ case 1 :
+ for (int i = 0; i < numEl; i++) {
+ int tag = method_1_2(el[i], maxDepth);
+ if (tag < 0) {
+ a[1]++;
+ if (tag < -1)
+ Msg::Info("Bad element : %d (with tag %d)", el[i]->getNum(), tag);
+ else
+ Msg::Info("Bad element : %d", el[i]->getNum());
+ }
+ else if (tag > 0) {
+ el[i]->setVisibility(0);
+ if (tag > 1)
+ Msg::Info("Good element : %d (with tag %d)", el[i]->getNum(), tag);
+ }
+ else {
+ a[0]++;
+ Msg::Info("Element %d may be bad", el[i]->getNum());
+ }
+ }
+ return a;
+
+ case 2 :
+ std::vector<int> tag(numEl, UNDEF_JAC_TAG);
+ method_2_2(el, tag, maxDepth);
+
+ Msg::Info(" ");
+ Msg::Info("Bad elements :");
+ for (unsigned int i = 0; i < tag.size(); i++) {
+ if (tag[i] < 0) {
+ if (tag[i] < -1)
+ Msg::Info("%d (with tag %d)", el[i]->getNum(), tag[i]);
+ else
+ Msg::Info("%d", el[i]->getNum());
+ a[1]++;
+ }
+ }
+
+ Msg::Info(" ");
+ Msg::Info("Uncertain elements :");
+ for (unsigned int i = 0; i < tag.size(); i++) {
+ if (tag[i] == 0) {
+ Msg::Info("%d", el[i]->getNum());
+ a[0]++;
+ }
+ }
+
+ Msg::Info(" ");
+ return a;
+ }
}
-int GMSH_AnalyseCurvedMeshPlugin::checkJacobian(MElement *el, int depth)
+
+int GMSH_AnalyseCurvedMeshPlugin::method_1_1(MElement *el, int depth)
{
- int tag = method1(el, depth);
- if (tag < 0) {
- Msg::Info("Bad element : %d", el->getNum());
- }
- else if (tag > 0) {
- el->setVisibility(0);
- }
- else {
- Msg::Info("Element %d may be bad.", el->getNum());
- }
- return tag;
+ const polynomialBasis *fs = el->getFunctionSpace(-1);
+ if (!fs) {
+ Msg::Error("Function space not implemented for type of element %d", el->getNum());
+ return 0;
+ }
+ const JacobianBasis *jfs = el->getJacobianFuncSpace(-1);
+ if (!jfs) {
+ Msg::Error("Jacobian function space not implemented for type of element %d", el->getNum());
+ return 0;
+ }
+
+ int numSamplingPt = jfs->points.size1();
+ int dim = jfs->points.size2();
+ fullVector<double> jacobian(numSamplingPt);
+
+ for (int i = 0; i < numSamplingPt; i++) {
+ double gsf[256][3];
+ switch (dim) {
+ case 1 :
+ fs->df(jfs->points(i,0),0,0, gsf);
+ break;
+ case 2 :
+ fs->df(jfs->points(i,0),jfs->points(i,1),0, gsf);
+ break;
+ case 3 :
+ fs->df(jfs->points(i,0),jfs->points(i,1),jfs->points(i,2), gsf);
+ break;
+ default :
+ Msg::Error("Can't get the gradient for %dD elements.", dim);
+ return false;
+ }
+ double jac[3][3];
+ jacobian(i) = getJacobian(gsf, jac, el);
+ }
+
+ for (int i = 0; i < jacobian.size(); i++) {
+ if (jacobian(i) <= 0.) return -1;
+ }
+
+
+ fullVector<double> jacBez(jacobian.size());
+ jfs->matrixLag2Bez.mult(jacobian, jacBez);
+
+ bool allPtPositive = true;
+ for (int i = 0; i < jacBez.size(); i++) {
+ if (jacBez(i) <= 0.) allPtPositive = false;
+ }
+ if (allPtPositive) return 1;
+
+
+ if (depth <= 1) {
+ return 0;
+ }
+ else{
+ int tag = division(jfs, jacBez, depth-1);
+ if (tag < 0)
+ return tag - 1;
+ if (tag > 0)
+ return tag + 1;
+ return tag;
+ }
}
-int GMSH_AnalyseCurvedMeshPlugin::method1(MElement *el, int depth)
+
+int GMSH_AnalyseCurvedMeshPlugin::method_1_2(MElement *el, int depth)
{
- const polynomialBasis *fs = el->getFunctionSpace(-1);
- if (!fs) {
- Msg::Error("Function space not implemented for type of element %d", el->getNum());
- return 0;
- }
- const JacobianBasis *jfs = el->getJacobianFuncSpace(-1);
- if (!jfs) {
- Msg::Error("Jacobian function space not implemented for type of element %d", el->getNum());
- return 0;
- }
- int numSamplingPt = jfs->points.size1(), dim = jfs->points.size2();
- fullVector<double> jacobian(numSamplingPt);
-
- for (int i = 0; i < numSamplingPt; i++) {
- double gsf[256][3];
- switch (dim) {
- case 1:
- fs->df(jfs->points(i,0),0,0, gsf);
- break;
- case 2:
- fs->df(jfs->points(i,0),jfs->points(i,1),0, gsf);
- break;
- case 3:
- fs->df(jfs->points(i,0),jfs->points(i,1),jfs->points(i,2), gsf);
- break;
- default:
- Msg::Error("Can't get the gradient for %dD elements.", dim);
- return false;
- }
- double jac[3][3];
- jacobian(i) = getJacobian(gsf, jac, el);
- }
-
- for (int i = 0; i < jacobian.size(); i++) {
- if (jacobian(i) < 0.) return -1;
- }
-
- fullVector<double> jacBez(jacobian.size());
- jfs->matrixLag2Bez.mult(jacobian, jacBez);
-
- bool allPtPositive = true;
- for (int i = 0; i < jacBez.size(); i++) {
- if (jacBez(i) < 0.) allPtPositive = false;
- }
- if (allPtPositive) return 1;
-
- if (depth <= 0) {
- return 0;
- }
- //else{return division(el, depth-1);}
-
- return 0;
-}
\ No newline at end of file
+ const polynomialBasis *fs = el->getFunctionSpace(-1);
+ if (!fs) {
+ Msg::Error("Function space not implemented for type of element %d", el->getNum());
+ return 0;
+ }
+ const JacobianBasis *jfs = el->getJacobianFuncSpace(-1);
+ if (!jfs) {
+ Msg::Error("Jacobian function space not implemented for type of element %d", el->getNum());
+ return 0;
+ }
+
+ int numSamplingPt = jfs->points.size1();
+ int dim = jfs->points.size2();
+ fullVector<double> jacobian(numSamplingPt);
+
+
+ setJacobian(el, jfs, jacobian);
+
+ //{
+ // Msg::Info("Printing vector jac[%d]", el->getNum());
+ // Msg::Info(" ");
+ // for(int I = 0; I < jacobian.size(); I++){
+ // Msg::Info("%lf ", jacobian(I));
+ // }
+ // Msg::Info(" ");
+ //}
+
+ for (int i = 0; i < jacobian.size(); i++) {
+ if (jacobian(i) <= 0.) return -1;
+ }
+
+
+ fullVector<double> jacBez(jacobian.size());
+ jfs->matrixLag2Bez.mult(jacobian, jacBez);
+
+ bool allPtPositive = true;
+ for (int i = 0; i < jacBez.size(); i++) {
+ if (jacBez(i) <= 0.) allPtPositive = false;
+ }
+ if (allPtPositive) return 1;
+
+
+ if (depth <= 1) {
+ return 0;
+ }
+ else{
+ int tag = division(jfs, jacBez, depth-1);
+ if (tag < 0)
+ return tag - 1;
+ if (tag > 0)
+ return tag + 1;
+ return tag;
+ }
+}
+
+
+void GMSH_AnalyseCurvedMeshPlugin::method_2_2
+ (MElement *const *el, std::vector<int> &tags, int depth)
+{
+ const polynomialBasis *fs = el[0]->getFunctionSpace(-1);
+ if (!fs) {
+ Msg::Error("Function space not implemented for type of element %d", el[0]->getNum());
+ return;
+ }
+ const JacobianBasis *jfs = el[0]->getJacobianFuncSpace(-1);
+ if (!jfs) {
+ Msg::Error("Jacobian function space not implemented for type of element %d", el[0]->getNum());
+ return;
+ }
+
+ int numSamplingPt = jfs->points.size1();
+ int dim = jfs->points.size2();
+ int numEl = tags.size();
+
+ fullMatrix<double> jacobian(numSamplingPt, numEl);
+ setJacobian(el, jfs, jacobian);
+
+
+ int numBad = 0;
+
+ for (int i = 0; i < numEl; i++) {
+ bool isBad = false;
+ for (int j = 0; j < jacobian.size1(); j++) {
+ if (jacobian(j,i) <= 0.) isBad = true;
+ }
+ if (isBad) {
+ tags[i] = -1;
+ numBad++;
+ }
+ }
+
+
+ fullMatrix<double> jacBez;
+ std::vector<int> correspondingEl;
+
+ double critere = .2; // à définir suivant le temps de chaque opération
+ if ((double) numBad / (double) numEl < critere) {// il vaut mieux calculer les points de controle de chaque élément
+ jacBez.resize(numSamplingPt, numEl);
+ jfs->matrixLag2Bez.mult(jacobian, jacBez);
+ correspondingEl.resize(numEl);
+ for (int i = 0; i < numEl; i++) correspondingEl[i] = i;
+ }
+ else {// il vaut mieux ne calculer que les points de controle des éléments encore incertains
+ fullMatrix<double> newJac(numSamplingPt, numEl-numBad);
+ fullMatrix<double> prox1(numSamplingPt,1);
+ fullMatrix<double> prox2(numSamplingPt,1);
+ int index = 0;
+ correspondingEl.resize(numEl - numBad);
+ for (int i = 0; i < numEl; i++) {
+ if (tags[i] == UNDEF_JAC_TAG) {
+ correspondingEl[index] = i;
+ prox1.setAsProxy(newJac,index,1);
+ prox2.setAsProxy(jacobian,i,1);
+ prox1.copy(prox2);
+ }
+ }
+ jacBez.resize(numSamplingPt, numEl-numBad);
+ jfs->matrixLag2Bez.mult(jacobian, jacBez);
+ }
+
+
+ int numGood = 0;
+
+ for (int i = 0; i < jacBez.size2(); i++) {
+ bool allPtPositive = true;
+ for (int j = 0; j < jacBez.size1(); j++) {
+ if (jacBez(j,i) <= 0.) allPtPositive = false;
+ }
+ if (allPtPositive) {
+ tags[correspondingEl[i]] = 1;
+ numGood++;
+ }
+ }
+
+
+ Msg::Warning("Not yet implemented for maxDepth > 1");
+ depth = 1;
+ if (depth <= 1) {
+ for (int i = 0; i < jacBez.size2(); i++)
+ if (tags[correspondingEl[i]] == UNDEF_JAC_TAG)
+ tags[correspondingEl[i]] = 0;
+ return;
+ }
+ /*else{
+ int tag = division(jfs, jacBez, depth-1);
+ if (tag < 0)
+ return tag - 1;
+ if (tag > 0)
+ return tag + 1;
+ return tag;
+ }*/
+
+}
+
+int GMSH_AnalyseCurvedMeshPlugin::division
+ (const JacobianBasis *jfs, const fullVector<double> &jacobian, int depth)
+{
+ if (jfs->divisor.size2() != jacobian.size()) {
+ Msg::Error("Wrong sizes in division : [%d,%d] * [%d]",
+ jfs->divisor.size1(), jfs->divisor.size2(), jacobian.size());
+ Msg::Info(" ");
+ return 0;
+ }
+
+ fullVector<double> newJacobian(jfs->divisor.size1());
+ jfs->divisor.mult(jacobian, newJacobian);
+
+ for (int i = 0; i < jfs->numDivisions; i++) {
+ for (int j = 0; j < jfs->numLagPts; j++)
+ if (newJacobian(i * jfs->points.size1() + j) <= 0.) return -1;
+ }
+
+ bool allPtPositive = true;
+ for (int i = 0; i < newJacobian.size(); i++) {
+ if (newJacobian(i) <= 0.) allPtPositive = false;
+ }
+ if (allPtPositive) return 1;
+
+
+ // We don't know if the jacobian is positif everywhere or not
+ // We subdivide if we still can do it (if depth > 0).
+ if (depth <= 0) {
+ return 0;
+ }
+ else{
+ fullVector<double> subJacobian;
+ std::vector<int> negTag, posTag;
+ bool zeroTag = false;
+
+ for (int i = 0; i < jfs->numDivisions; i++)
+ {
+ subJacobian.setAsProxy(newJacobian, i * jacobian.size(), jacobian.size());
+ int tag = division(jfs, subJacobian, depth-1);
+
+ if (tag < 0)
+ negTag.push_back(tag);
+ else if (tag > 0)
+ posTag.push_back(tag);
+ else
+ zeroTag = true;
+ }
+
+ if (negTag.size() > 0) return max(negTag) - 1;
+
+ if (zeroTag) return 0;
+
+ return max(posTag) - 1;
+ }
+}
+
+
+
+
+ /*{
+ Msg::Info("Printing matrix %s :", "nodesX");
+ int ni = nodesX.size1();
+ int nj = nodesX.size2();
+ for(int I = 0; I < ni; I++){
+ Msg::Info(" ");
+ for(int J = 0; J < nj; J++){
+ Msg::Info("%lf ", nodesX(I, J));
+ }
+ }
+ Msg::Info(" ");
+ }*/
\ No newline at end of file
diff --git a/Plugin/AnalyseCurvedMesh.h b/Plugin/AnalyseCurvedMesh.h
index cf74b0722bc734d5f45698b404d78dc92bf89ff3..985cdc73835e64ad70a78bea507c16f7e373958f 100644
--- a/Plugin/AnalyseCurvedMesh.h
+++ b/Plugin/AnalyseCurvedMesh.h
@@ -2,20 +2,13 @@
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
-//
-// Contributed by Amaury J.
-
-// Added code :
-// Numeric/polynomialBasis.cpp -> jacobianPolynomialBases::find(int tag), #include
-// Geo/MElement.h -> getJacobianFunctionSpace(int), #include
-// Geo/MTriangle .h/.cpp -> getJacobianFunctionSpace(int)
#ifndef _ANALYSECURVEDMESH_H_
#define _ANALYSECURVEDMESH_H_
-#include <string>
#include "Plugin.h"
#include "MElement.h"
+#include <vector>
extern "C"
{
@@ -33,10 +26,19 @@ class GMSH_AnalyseCurvedMeshPlugin : public GMSH_PostPlugin
}
std::string getHelp() const;
std::string getAuthor() const { return "Amaury Johnen"; }
+ int getNbOptions() const;
+ StringXNumber *getOption(int);
PView *execute(PView *);
bool isJacPositive(MElement *);
- int method1(MElement *, int depth);
- int checkJacobian(MElement *, int depth);
+ int method_1_1(MElement *, int depth);
+ int method_1_2(MElement *, int depth);
+ int method_1_3(MElement *, int depth);
+ void method_2_2(MElement *const *, std::vector<int> &tags, int depth);
+ void method_2_3(MElement *const *, std::vector<int> &tags, int depth);
+ //int checkJacobian(MElement *, int depth);
+ //int *checkJacobian2(MElement *const *, int numEl, int depth);
+ int *checkJacobian(MElement *const *, int numEl, int depth, int method);
+ int division(const JacobianBasis *, const fullVector<double> &, int depth);
};
#endif