From 4d2b0ea1747d2d09d3dd388ccb9d51ca21f55db1 Mon Sep 17 00:00:00 2001
From: Amaury Johnan <amjohnen@gmail.com>
Date: Wed, 4 Aug 2010 11:59:49 +0000
Subject: [PATCH] AnalyseCurvedMesh development (works for triangle and tetra
but is not yet optimized and is still incomplete, edit PluginManager.cpp to
activate the plugin)
---
Geo/MElement.h | 4 +
Geo/MTetrahedron.cpp | 34 +
Geo/MTetrahedron.h | 1 +
Geo/MTriangle.cpp | 39 +
Geo/MTriangle.h | 1 +
Numeric/CMakeLists.txt | 1 +
Numeric/JacobianBasis.cpp | 1679 ++++++++++++++++++++++++++++++++++
Numeric/JacobianBasis.h | 237 +++++
Numeric/fullMatrix.cpp | 4 +-
Numeric/fullMatrix.h | 2 +-
Plugin/AnalyseCurvedMesh.cpp | 332 ++++++-
Plugin/AnalyseCurvedMesh.h | 11 +-
12 files changed, 2293 insertions(+), 52 deletions(-)
create mode 100644 Numeric/JacobianBasis.cpp
create mode 100644 Numeric/JacobianBasis.h
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 3e9494ca8b..55e7011c03 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -15,6 +15,7 @@
#include "MEdge.h"
#include "MFace.h"
#include "polynomialBasis.h"
+#include "JacobianBasis.h"
#include "Gauss.h"
class GFace;
@@ -205,6 +206,9 @@ class MElement
// get the function space for the element
virtual const polynomialBasis* getFunctionSpace(int order=-1) const { return 0; }
+ // get the function space for the jacobian of the element
+ virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const { return 0; }
+
// return the interpolating nodal shape functions evaluated at point
// (u,v,w) in parametric coordinates (if order == -1, use the
// polynomial order of the element)
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index 2cdeb6704e..2612ac62ed 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -137,6 +137,40 @@ const polynomialBasis* MTetrahedron::getFunctionSpace(int o) const
return 0;
}
+const JacobianBasis* MTetrahedron::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_TET_4);
+ case 2: return JacobianBases::find(MSH_TET_10);
+ case 3: return JacobianBases::find(MSH_TET_20);
+ case 4: return JacobianBases::find(MSH_TET_34);
+ case 5: return JacobianBases::find(MSH_TET_52);
+ default: Msg::Error("Order %d tetrahedron function space not implemented", order);
+ }
+ }
+ else {
+ switch (order) {
+ case 1: return JacobianBases::find(MSH_TET_4);
+ case 2: return JacobianBases::find(MSH_TET_10);
+ case 3: return JacobianBases::find(MSH_TET_20);
+ case 4: return JacobianBases::find(MSH_TET_35);
+ case 5: return JacobianBases::find(MSH_TET_56);
+ case 6: return JacobianBases::find(MSH_TET_84);
+ case 7: return JacobianBases::find(MSH_TET_120);
+ case 8: return JacobianBases::find(MSH_TET_165);
+ case 9: return JacobianBases::find(MSH_TET_220);
+ case 10: return JacobianBases::find(MSH_TET_286);
+ default: Msg::Error("Order %d tetrahedron function space not implemented", order);
+ }
+ }
+ return 0;
+}
+
int MTetrahedron10::getNumEdgesRep(){ return 6 * CTX::instance()->mesh.numSubEdges; }
int MTetrahedronN::getNumEdgesRep(){ return 6 * CTX::instance()->mesh.numSubEdges; }
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index bd817daae5..e0877afd6f 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -133,6 +133,7 @@ class MTetrahedron : public MElement {
virtual double etaShapeMeasure();
void xyz2uvw(double xyz[3], double uvw[3]);
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/MTriangle.cpp b/Geo/MTriangle.cpp
index 78fc812fa1..90493793f5 100644
--- a/Geo/MTriangle.cpp
+++ b/Geo/MTriangle.cpp
@@ -104,6 +104,45 @@ const polynomialBasis* MTriangle::getFunctionSpace(int o) const
return 0;
}
+const JacobianBasis* MTriangle::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_TRI_3);
+ case 2: return JacobianBases::find(MSH_TRI_6);
+ case 3: return JacobianBases::find(MSH_TRI_9);
+ case 4: return JacobianBases::find(MSH_TRI_12);
+ case 5: return JacobianBases::find(MSH_TRI_15I);
+ case 6: return JacobianBases::find(MSH_TRI_18);
+ case 7: return JacobianBases::find(MSH_TRI_21I);
+ case 8: return JacobianBases::find(MSH_TRI_24);
+ case 9: return JacobianBases::find(MSH_TRI_27);
+ case 10: return JacobianBases::find(MSH_TRI_30);
+ default: Msg::Error("Order %d triangle incomplete function space not implemented", order);
+ }
+ }
+ else {
+ switch (order) {
+ case 1: return JacobianBases::find(MSH_TRI_3);
+ case 2: return JacobianBases::find(MSH_TRI_6);
+ case 3: return JacobianBases::find(MSH_TRI_10);
+ case 4: return JacobianBases::find(MSH_TRI_15);
+ case 5: return JacobianBases::find(MSH_TRI_21);
+ case 6: return JacobianBases::find(MSH_TRI_28);
+ case 7: return JacobianBases::find(MSH_TRI_36);
+ case 8: return JacobianBases::find(MSH_TRI_45);
+ case 9: return JacobianBases::find(MSH_TRI_55);
+ case 10: return JacobianBases::find(MSH_TRI_66);
+ default: Msg::Error("Order %d triangle function space not implemented", order);
+ }
+ }
+ return 0;
+}
+
int MTriangleN::getNumEdgesRep(){ return 3 * CTX::instance()->mesh.numSubEdges; }
int MTriangle6::getNumEdgesRep(){ return 3 * CTX::instance()->mesh.numSubEdges; }
diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h
index 091c583319..53faaf9e32 100644
--- a/Geo/MTriangle.h
+++ b/Geo/MTriangle.h
@@ -126,6 +126,7 @@ class MTriangle : public MElement {
MVertex *tmp = _v[1]; _v[1] = _v[2]; _v[2] = 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/Numeric/CMakeLists.txt b/Numeric/CMakeLists.txt
index 3a9b71fe38..edf76e892e 100644
--- a/Numeric/CMakeLists.txt
+++ b/Numeric/CMakeLists.txt
@@ -7,6 +7,7 @@ set(SRC
Numeric.cpp
fullMatrix.cpp
polynomialBasis.cpp
+ JacobianBasis.cpp
GaussQuadratureLin.cpp
GaussQuadratureTri.cpp
GaussQuadratureQuad.cpp
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
new file mode 100644
index 0000000000..4ff2f8b933
--- /dev/null
+++ b/Numeric/JacobianBasis.cpp
@@ -0,0 +1,1679 @@
+// 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
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
new file mode 100644
index 0000000000..885cf3b1a1
--- /dev/null
+++ b/Numeric/JacobianBasis.h
@@ -0,0 +1,237 @@
+// 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>.
+
+#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;
+};
+
+class JacobianBases
+{
+ private:
+ static std::map<int, JacobianBasis> fs;
+ public:
+ 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 c5924109c7..faa67bd3ec 100644
--- a/Numeric/fullMatrix.cpp
+++ b/Numeric/fullMatrix.cpp
@@ -110,7 +110,7 @@ void fullMatrix<std::complex<double> >::gemm(const fullMatrix<std::complex<doubl
}
template<>
-void fullMatrix<double>::mult(const fullVector<double> &x, fullVector<double> &y)
+void fullMatrix<double>::mult(const fullVector<double> &x, fullVector<double> &y) const
{
int M = _r, N = _c, LDA = _r, INCX = 1, INCY = 1;
double alpha = 1., beta = 0.;
@@ -120,7 +120,7 @@ void fullMatrix<double>::mult(const fullVector<double> &x, fullVector<double> &y
template<>
void fullMatrix<std::complex<double> >::mult(const fullVector<std::complex<double> > &x,
- fullVector<std::complex<double> > &y)
+ fullVector<std::complex<double> > &y) const
{
int M = _r, N = _c, LDA = _r, INCX = 1, INCY = 1;
std::complex<double> alpha = 1., beta = 0.;
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index cb43fa928d..aedad7f1e3 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -321,7 +321,7 @@ class fullMatrix
for(int j = 0; j < size2(); j++)
(*this)(i, j) += a*m(i, j);
}
- void mult(const fullVector<scalar> &x, fullVector<scalar> &y)
+ void mult(const fullVector<scalar> &x, fullVector<scalar> &y) const
#if !defined(HAVE_BLAS)
{
y.scale(0.);
diff --git a/Plugin/AnalyseCurvedMesh.cpp b/Plugin/AnalyseCurvedMesh.cpp
index 3230647024..0f0d3ba6ec 100644
--- a/Plugin/AnalyseCurvedMesh.cpp
+++ b/Plugin/AnalyseCurvedMesh.cpp
@@ -1,49 +1,285 @@
-// 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>.
-//
-// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
+// 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>.
+//
+// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
-#include <stdlib.h>
-#include "Gmsh.h"
-#include "MTriangle.h"
-#include "GmshConfig.h"
-#include "GModel.h"
-#include "polynomialBasis.h"
-#include "AnalyseCurvedMesh.h"
-
-
-extern "C"
-{
- GMSH_Plugin *GMSH_RegisterAnalyseCurvedMeshPlugin()
- {
- return new GMSH_AnalyseCurvedMeshPlugin();
- }
-}
-
-std::string GMSH_AnalyseCurvedMeshPlugin::getHelp() const
-{
- return "Plugin(AnalyseCurvedMesh) computes AnalyseCurvedMeshs to boundaries in "
- "a mesh.\n\n"
- "Plugin(AnalyseCurvedMesh) creates a bunch of files.";
-}
-
-PView *GMSH_AnalyseCurvedMeshPlugin::execute(PView *v)
-{
- Msg::Info("Hello !");
-
- GModel *m = GModel::current();
-
- for(GModel::fiter it = m->firstFace(); it != m->lastFace(); it++) {
- GFace *f = *it;
- Msg::Info("Je suis dans la surface %d", f->tag());
-
- for(int i = 0; i < (int) f->triangles.size(); i++) {
- MTriangle *t = f->triangles[i];
- const polynomialBasis *p = t->getFunctionSpace();
- Msg::Info("Taille de la matrice P : %dx%d", p->monomials.size1(), p->monomials.size2());
- }
- }
- return 0;
-}
+#include "GmshDefines.h"
+#include <stdlib.h>
+#include "Gmsh.h"
+#include "MTriangle.h"
+#include "GmshConfig.h"
+#include "GModel.h"
+#include "polynomialBasis.h"
+#include "JacobianBasis.h"
+#include "AnalyseCurvedMesh.h"
+#include "fullMatrix.h"
+
+extern "C"
+{
+ GMSH_Plugin *GMSH_RegisterAnalyseCurvedMeshPlugin()
+ {
+ return new GMSH_AnalyseCurvedMeshPlugin();
+ }
+}
+
+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)";
+}
+
+static double computeDeterminant(MElement *el, double jac[3][3])
+{
+ switch (el->getDim()) {
+ case 0:
+ return 1.0;
+ case 1:
+ return jac[0][0];
+ case 2:
+ return jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0];
+ case 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;
+ }
+}
+
+double getJacobian(double gsf[][3], double jac[3][3], MElement *el)
+{
+ jac[0][0] = jac[0][1] = jac[0][2] = 0.;
+ jac[1][0] = jac[1][1] = jac[1][2] = 0.;
+ jac[2][0] = jac[2][1] = jac[2][2] = 0.;
+
+ for (int i = 0; i < el->getNumVertices(); i++) {
+ const MVertex *v = el->getVertex(i);
+ double *gg = gsf[i];
+ for (int j = 0; j < 3; j++) {
+ jac[j][0] += v->x() * gg[j];
+ jac[j][1] += v->y() * gg[j];
+ jac[j][2] += v->z() * gg[j];
+ }
+ }
+
+ return computeDeterminant(el, jac);
+}
+
+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.");
+
+ return 0;
+}
+
+bool GMSH_AnalyseCurvedMeshPlugin::isJacPositive(MElement *el)
+{
+ 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 GMSH_AnalyseCurvedMeshPlugin::checkJacobian(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;
+}
+
+int GMSH_AnalyseCurvedMeshPlugin::method1(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
diff --git a/Plugin/AnalyseCurvedMesh.h b/Plugin/AnalyseCurvedMesh.h
index a67da0b067..cf74b0722b 100644
--- a/Plugin/AnalyseCurvedMesh.h
+++ b/Plugin/AnalyseCurvedMesh.h
@@ -3,13 +3,19 @@
// 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>.
+// 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"
extern "C"
{
@@ -28,6 +34,9 @@ class GMSH_AnalyseCurvedMeshPlugin : public GMSH_PostPlugin
std::string getHelp() const;
std::string getAuthor() const { return "Amaury Johnen"; }
PView *execute(PView *);
+ bool isJacPositive(MElement *);
+ int method1(MElement *, int depth);
+ int checkJacobian(MElement *, int depth);
};
#endif
--
GitLab