diff --git a/Common/Gmsh.cpp b/Common/Gmsh.cpp
index 4bb15f44193a3fbe09ede3e2c9f9289bb7182f11..e97c8396916ce655f566eb46b4afdb3d8b1aee4f 100644
--- a/Common/Gmsh.cpp
+++ b/Common/Gmsh.cpp
@@ -244,89 +244,6 @@ int GmshBatch()
solver_batch_cb(0, (void*)CTX::instance()->launchSolverAtStartup);
#endif
- if (false) {
- // 11/06/13 : This will be removed later !
- std::vector<int> compare;
- std::vector<int> test;
- for (int i = 1; i < MSH_NUM_TYPE+1; ++i) {
- if (i == 76 || i == 77 || i == 78)
- continue;
-
- const char **name = new const char*[1];
- MElement::getInfoMSH(i, name);
-
- if (i == 67 || i == 68 || i == 70) {
- Msg::Info("ignoring unknown %d", i);
- continue;
- }
- if (MElement::ParentTypeFromTag(i) == TYPE_POLYG) {
- Msg::Info("ignoring polygone %d", i);
- continue;
- }
- if (MElement::ParentTypeFromTag(i) == TYPE_POLYH) {
- Msg::Info("ignoring polyhedre %d", i);
- continue;
- }
- if (MElement::ParentTypeFromTag(i) == TYPE_XFEM) {
- Msg::Info("ignoring xfem %d", i);
- continue;
- }
-
- if (MElement::ParentTypeFromTag(i) == TYPE_PRI && MElement::OrderFromTag(i) > 2) {
- Msg::Warning("%d:ignoring %s (different node order)", i, *name);
- test.push_back(i);
- }
- else if (MElement::ParentTypeFromTag(i) == TYPE_PYR && MElement::SerendipityFromTag(i) > 1) {
- Msg::Warning("%d:ignoring Serendipity %s (bug in Bergot implementation & new Basis function not implemented)", i, *name);
- }
- else if (MElement::DimensionFromTag(i) == 3 && MElement::SerendipityFromTag(i) > 1) {
- Msg::Warning("%d:ignoring 3D Serendipity : %s (different definition)", i, *name);
- test.push_back(i);
- }
- else {
- const nodalBasis *fs = BasisFactory::getNodalBasis(i);
- if (fs && !fs->compareNewAlgoPointsWithOld()) {
- compare.push_back(i);
- }
- }
- }
- Msg::Info(" ");
- for (int parent = 1; parent < 12; ++parent) {
- for (unsigned int i = 0; i < compare.size(); ++i) {
- if (MElement::ParentTypeFromTag(compare[i]) == parent) {
- if (parent == TYPE_PYR) {
- const char **name = new const char*[1];
- MElement::getInfoMSH(compare[i], name);
- Msg::Warning("%d, %s: ignoring, old implementation much better and will be kept", compare[i], *name);
- }
- else {
- const nodalBasis *fs = BasisFactory::getNodalBasis(compare[i]);
- fs->compareNewAlgoBaseFunctionsWithOld();
- }
- }
- }
- }
- Msg::Info(" ");
- Msg::Info("---------------------------");
- for (unsigned int i = 0; i < test.size(); ++i) {
- const char **name = new const char*[1];
- MElement::getInfoMSH(test[i], name);
- //Msg::Info("%d, %s:", test[i], *name);
- const nodalBasis *fs = BasisFactory::getNodalBasis(test[i]);
- }
- Msg::Info("---------------------------");
- for (int parent = 1; parent < 12; ++parent) {
- for (unsigned int i = 0; i < test.size(); ++i) {
- if (MElement::ParentTypeFromTag(test[i]) == parent) {
- const nodalBasis *fs = BasisFactory::getNodalBasis(test[i]);
- /*const char **name = new const char*[1];
- MElement::getInfoMSH(test[i], name);
- Msg::Warning("%d, %s: untested !", test[i], *name);*/
- fs->testNewAlgoBaseFunctions();
- }
- }
- }
- }
time_t now;
time(&now);
diff --git a/Numeric/nodalBasis.cpp b/Numeric/nodalBasis.cpp
index 2e1a2d0db543907ab0f20d3e95cefa428746cafe..da3de9954cd8932368e9b0c56748c09f753d3407 100644
--- a/Numeric/nodalBasis.cpp
+++ b/Numeric/nodalBasis.cpp
@@ -10,155 +10,6 @@
//#include "pointsGenerators.h"
-int nodalBasis::compareNewAlgoPointsWithOld() const
-{
- const char **name = new const char*[1];
- MElement::getInfoMSH(type, name);
- if (points_newAlgo.size1() == 0) {
- Msg::Error("%d: pas de points (%d, %d, %d) %s", type, parentType, serendip, order, *name);
- return 1;
- }
- if (points_newAlgo.size1() != points.size1()) {
- Msg::Error("%d: pas meme taille (%d, %d, %d) %s", type, parentType, serendip, order, *name);
- Msg::Error(" |--> points[%d, %d] vs newPoints[%d, %d]", points.size1(), points.size2(), points_newAlgo.size1(), points_newAlgo.size2());
- return 2;
- }
- for (int i = 0; i < points.size1(); ++i) {
- for (int j = 0; j < points.size2(); ++j) {
- //Msg::Info("(i, j) = (%d, %d)", i, j);
- //Msg::Info(" ");
- if (std::abs(points(i, j) - points_newAlgo(i, j)) > 1e-15) {
- Msg::Error("%d: correspond pas (%d, %d, %d) %s", type, parentType, serendip, order, *name);
- for (int i = 0; i < points.size1(); ++i) {
- for (int j = 0; j < points.size2(); ++j) {
- if(std::abs(points(i, j) - points_newAlgo(i, j)) <= 1e-15)
- Msg::Info("(%d,%d) : points %f / %f newPoints (mon %d)", i, j, points(i, j), points_newAlgo(i, j), monomials_newAlgo(i, j));
- else
- Msg::Error("(%d,%d) : points %f / %f newPoints (mon %d)", i, j, points(i, j), points_newAlgo(i, j), monomials_newAlgo(i, j));
- }
- Msg::Info(" ");
- }
- return 3;
- }
- }
- }
- Msg::Info("%d: ok (%d, %d, %d) %s", type, parentType, serendip, order, *name);
- return 0;
-}
-
-int nodalBasis::compareNewAlgoBaseFunctionsWithOld() const
-{
- const char **name = new const char*[1];
- MElement::getInfoMSH(type, name);
- int ndof = points.size1();
- int P[3];
- double oldVal[ndof], newVal[ndof];
- for (int i = 0; i < 10; ++i) {
- if (i == 0) {
- P[0] = 0;
- P[1] = 0;
- P[2] = 0;
- }
- else if (i == 1) {
- P[0] = 0;
- P[1] = 0;
- P[2] = 1000;
- }
- else if (i == 2) {
- P[0] = 1000;
- P[1] = 0;
- P[2] = 0;
- }
- else {
- P[0] = std::rand() % 1001;
- P[1] = std::rand() % (1001 - P[0]);
- P[2] = std::rand() % (1001 - P[0] - P[1]);
- }
-
- f(P[0] / 1000., P[1] / 1000., P[2] / 1000., oldVal);
- fnew(P[0] / 1000., P[1] / 1000., P[2] / 1000., newVal);
-
- double sumError = .0;
- double sumOne[2];
- sumOne[0] = .0;
- sumOne[1] = .0;
- for (int j = 0; j < ndof; ++j) {
- double diff = std::abs(oldVal[j] - newVal[j]);
- //double mean = std::sqrt(.5 * (oldVal[j]*oldVal[j] + newVal[j]*newVal[j]));
- //double error = diff/(mean + 1e-15);
- double error = diff;
- sumError += error*error;
- sumOne[0] += oldVal[j];
- sumOne[1] += newVal[j];
- }
- sumError = std::sqrt(sumError/ndof);
- /*if (sumError > 1e-5) {
- Msg::Error("(%.2f, %.2f, %.2f) -> fold=%.2e / fnew=%.2e", P[0] / 1000., P[1] / 1000., P[2] / 1000., sumOne[0], sumOne[1]);
- //for (int j = 0; j < ndof; ++j) {
- // Msg::Info("old %.3e vs %.3e new", oldVal[j], newVal[j]);
- //}
- }*/
- if (sumError > 1e-5) {
- Msg::Error("%d, %s: f(%.2f, %.2f, %.2f) bad precision diff %.2e", type, *name, P[0] / 1000., P[1] / 1000., P[2] / 1000., sumError);
- return 1;
- }
- else if (sumError > 1e-13) {
- Msg::Warning("%d, %s: f(%.2f, %.2f, %.2f) bad precision diff %.2e", type, *name, P[0] / 1000., P[1] / 1000., P[2] / 1000., sumError);
- return 1;
- }
- }
- return 0;
-}
-
-int nodalBasis::testNewAlgoBaseFunctions() const
-{
- const char **name = new const char*[1];
- MElement::getInfoMSH(type, name);
- int ndof = points_newAlgo.size1();
- int P[3];
- bool noproblem = true;
- for (int i = 0; i < 10; ++i) {
- if (i == 0) {
- P[0] = 0;
- P[1] = 0;
- P[2] = 0;
- }
- else if (i == 1) {
- P[0] = 0;
- P[1] = 0;
- P[2] = 1000;
- }
- else if (i == 2) {
- P[0] = 1000;
- P[1] = 0;
- P[2] = 0;
- }
- else {
- P[0] = std::rand() % 1001;
- P[1] = std::rand() % (1001 - P[0]);
- P[2] = std::rand() % (1001 - P[0] - P[1]);
- }
-
- double newVal[ndof];
- fnew(P[0] / 1000., P[1] / 1000., P[2] / 1000., newVal);
-
- double sumOne = 0;
- for (int j = 0; j < ndof; ++j) {
- sumOne += newVal[j];
- }
- if (sumOne > 1. + 1e-12 || sumOne < 1. - 1e-12) {
- noproblem = false;
- Msg::Warning("%d, %s: bad precision (%.2f, %.2f, %.2f) -> sum = %.2e (1 + %.2e)", type, *name, P[0] / 1000., P[1] / 1000., P[2] / 1000., sumOne, sumOne-1.);
- //for (int j = 0; j < ndof; ++j) {
- // Msg::Info(" %d : %.3e", j, newVal[j]);
- //}
- return 1;
- }
- }
- if (noproblem) Msg::Info("%d, %s: no problem", type, *name);
- return 0;
-}
-
static int nbdoftriangle(int order) { return (order + 1) * (order + 2) / 2; }
//KH : caveat : node coordinates are not yet coherent with node numbering associated
diff --git a/Numeric/nodalBasis.h b/Numeric/nodalBasis.h
index d9d48e911b79a7f8d0a7b5bf45d311cccb543da5..ead6c1ef4f20f4b986b30b13bd68c6a392c4b179 100644
--- a/Numeric/nodalBasis.h
+++ b/Numeric/nodalBasis.h
@@ -32,10 +32,6 @@ class nodalBasis {
virtual void ddf(double u, double v, double w, double grads[][3][3]) const {Msg::Fatal("Not implemented");};
virtual void dddf(double u, double v, double w, double grads[][3][3][3]) const {Msg::Fatal("Not implemented");};
- int compareNewAlgoPointsWithOld() const;
- int compareNewAlgoBaseFunctionsWithOld() const;
- int testNewAlgoBaseFunctions() const;
-
// closures is the list of the nodes of each face, in the order of
// the polynomialBasis of the face; fullClosures is mapping of the
// nodes of the element that rotates the element so that the
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 5a7b626b9b19ced3fbdd8f3121678f92b4a081c5..24826d3845200b1b08758b4a601bc1cd7464684f 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -423,18 +423,6 @@ static fullMatrix<double> generateLagrangeMonomialCoefficients
fullMatrix<double> coefficient(ndofs, ndofs);
Vandermonde.invert(coefficient);
- fullMatrix<double> unity(ndofs, ndofs);
- Vandermonde.mult(coefficient, unity);
- double max = .0;
- for (int i = 0; i < ndofs; i++) {
- for (int j = 0; j < ndofs; j++) {
- if (i == j) unity(i, j) -= 1.;
- //Msg::Info(" unity(%d, %d) = %.3e", i, j, unity(i, j));
- max = std::max(max, std::abs(unity(i, j)));
- }
- }
- //if (max > 1e-10) Msg::Info(" max unity = %.3e", max);
-
return coefficient;
}
diff --git a/Numeric/pyramidalBasis.cpp b/Numeric/pyramidalBasis.cpp
index 1ee2fd4a0a0b9111e89c4091a0c94e53e17cbeac..2a3efeea5ea6c2dfcfae95cd48b58dcdaf310b5d 100644
--- a/Numeric/pyramidalBasis.cpp
+++ b/Numeric/pyramidalBasis.cpp
@@ -52,18 +52,6 @@ static fullMatrix<double> generateLagrangeMonomialCoefficientsPyr
fullMatrix<double> coefficient(ndofs, ndofs);
Vandermonde.invert(coefficient);
- fullMatrix<double> unity(ndofs, ndofs);
- Vandermonde.mult(coefficient, unity);
- double max = .0;
- for (int i = 0; i < ndofs; i++) {
- for (int j = 0; j < ndofs; j++) {
- if (i == j) unity(i, j) -= 1.;
- //Msg::Info(" unity(%d, %d) = %.3e", i, j, unity(i, j));
- max = std::max(max, std::abs(unity(i, j)));
- }
- }
- //if (max > 1e-14) Msg::Info("mon max unity = %.3e", max);
-
return coefficient;
}
@@ -88,19 +76,6 @@ pyramidalBasis::pyramidalBasis(int tag) : nodalBasis(tag)
delete[] fval;
- fullMatrix<double> unity(num_points, num_points);
- VDM.mult(VDMinv, unity);
- double max = .0;
- for (int i = 0; i < num_points; i++) {
- for (int j = 0; j < num_points; j++) {
- if (i == j) unity(i, j) -= 1.;
- //Msg::Info(" unity(%d, %d) = %.3e", i, j, unity(i, j));
- max = std::max(max, std::abs(unity(i, j)));
- }
- }
- //if (max > 1e-14) Msg::Info("leg max unity = %.3e", max);
-
-
// TEST NEW ALGO POINTS / MONOMIAL
monomials_newAlgo = gmshGenerateMonomialsPyramid(order, serendip);