diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 4b7ac5cb975f38dbac61905dcd8d3c5ea7f8ffaa..14a0689ed96fffd95508dc32b17391f924c29541 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -291,7 +291,7 @@ double MElement::getJacobian(const fullMatrix<double> &gsf, double jac[3][3])
for (int i = 0; i < getNumShapeFunctions(); i++) {
const MVertex *v = getShapeFunctionNode(i);
- for (int j = 0; j < 3; j++) {
+ for (int j = 0; j < gsf.size2(); j++) {
jac[j][0] += v->x() * gsf(i, j);
jac[j][1] += v->y() * gsf(i, j);
jac[j][2] += v->z() * gsf(i, j);
@@ -1216,70 +1216,6 @@ MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
}
}
-const fullMatrix<double> &MElement::
-getGradShapeFunctionsAtIntegrationPoints(int integrationOrder, int functionSpaceOrder)
-{
- static std::map <std::pair<int,int>, fullMatrix<double> > DF;
- const polynomialBasis *fs = getFunctionSpace(functionSpaceOrder);
- fullMatrix<double> &mat = DF[std::make_pair(fs->type, integrationOrder)];
- if (mat.size1() != 0) return mat;
- int npts;
- IntPt *pts;
- getIntegrationPoints(integrationOrder, &npts, &pts);
- mat.resize(fs->points.size1(), npts*3);
- double df[512][3];
- for (int i = 0; i < npts; i++) {
- fs->df(pts[i].pt[0], pts[i].pt[1], pts[i].pt[2], df);
- for (int j = 0; j < fs->points.size1(); j++) {
- mat(j, i*3+0) = df[j][0];
- mat(j, i*3+1) = df[j][1];
- mat(j, i*3+2) = df[j][2];
- }
- }
- return mat;
-}
-
-const fullMatrix<double> &MElement::
-getShapeFunctionsAtIntegrationPoints (int integrationOrder, int functionSpaceOrder)
-{
- static std::map <std::pair<int,int>, fullMatrix<double> > F;
- const polynomialBasis *fs = getFunctionSpace (functionSpaceOrder);
- fullMatrix<double> &mat = F[std::make_pair(fs->type, integrationOrder)];
- if (mat.size1() != 0) return mat;
- int npts;
- IntPt *pts;
- getIntegrationPoints(integrationOrder, &npts, &pts);
- mat.resize(fs->points.size1(), npts*3);
- double f[512];
- for (int i = 0; i < npts; i++) {
- fs->f(pts[i].pt[0], pts[i].pt[1], pts[i].pt[2], f);
- for (int j = 0; j < fs->points.size1(); j++) {
- mat(j, i) = f[j];
- }
- }
- return mat;
-}
-
-const fullMatrix<double> &MElement::getGradShapeFunctionsAtNodes(int functionSpaceOrder)
-{
- static std::map <int, fullMatrix<double> > DF;
- const polynomialBasis *fs = getFunctionSpace(functionSpaceOrder);
- fullMatrix<double> &mat = DF[fs->type];
- if (mat.size1()!=0) return mat;
- const fullMatrix<double> &points = fs->points;
- mat.resize (points.size1(), points.size1()*3);
- double df[512][3];
- for (int i = 0; i < points.size1(); i++) {
- fs->df (points(i,0), points(i,1), points(i,2), df);
- for (int j = 0; j < points.size1(); j++) {
- mat(j, i*3+0) = df[j][0];
- mat(j, i*3+1) = df[j][1];
- mat(j, i*3+2) = df[j][2];
- }
- }
- return mat;
-}
-
void MElement::xyzTouvw(fullMatrix<double> *xu)
{
double _xyz[3] = {(*xu)(0,0),(*xu)(0,1),(*xu)(0,2)}, _uvw[3];
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 534f22c1a4096d96351008e6a5e1e1a589b061dd..ae035be37199200eb2add8d566b3bd265446d5ed 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -236,12 +236,6 @@ class MElement
int order=-1);
virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3],
int order=-1);
- const fullMatrix<double> &getGradShapeFunctionsAtIntegrationPoints
- (int integrationOrder, int functionSpaceOrder=-1);
- const fullMatrix<double> &getShapeFunctionsAtIntegrationPoints
- (int integrationOrder, int functionSpaceOrder=-1);
- const fullMatrix<double> &getGradShapeFunctionsAtNodes (int functionSpaceOrder=-1);
-
// return the Jacobian of the element evaluated at point (u,v,w) in
// parametric coordinates
double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]);
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 6d7eac566beaeee441174ee64b5049d2717957d0..18b2caf491d553d135483d7afdc16629baa8c3fc 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -191,43 +191,6 @@ static fullMatrix<double> generatePascalSerendipityTriangle(int order)
return monomials;
}
-const fullMatrix<double> &polynomialBasis::
-getGradientAtFaceIntegrationPoints(int integrationOrder, int closureId) const
-{
- std::vector<fullMatrix<double> > &dfAtFace = _dfAtFace[integrationOrder];
- if (dfAtFace.empty()) {
- dfAtFace.resize(closures.size());
- int nbPsi = points.size1();
- for (unsigned int iClosure = 0; iClosure < closures.size(); iClosure++) {
- fullMatrix<double> integrationFace, weight;
- const polynomialBasis *fsFace = polynomialBases::find(closures[iClosure].type);
- gaussIntegration::get(fsFace->parentType, integrationOrder, integrationFace, weight);
- fullMatrix<double> integration(integrationFace.size1(), 3);
- double f[1256];
- for (int i = 0; i < integrationFace.size1(); i++){
- fsFace->f(integrationFace(i,0),integrationFace(i,1),integrationFace(i,2),f);
- for(unsigned int j = 0; j < closures[iClosure].size(); j++) {
- int jNod = closures[iClosure][j];
- for(int k = 0; k < points.size2(); k++)
- integration(i,k) += f[j] * points(jNod,k);
- }
- }
- fullMatrix<double> &v = dfAtFace[iClosure];
- v.resize(nbPsi, integration.size1()*3);
- double g[1256][3];
- for (int xi = 0; xi < integration.size1(); xi++) {
- df(integration(xi, 0), integration(xi, 1), integration(xi, 2), g);
- for (int j = 0; j < nbPsi; j++) {
- v(j, xi*3) = g[j][0];
- v(j, xi*3+1) = g[j][1];
- v(j, xi*3+2) = g[j][2];
- }
- }
- }
- }
- return dfAtFace[closureId];
-}
-
// generate all monomials xi^m * eta^n with n and m <= order
static fullMatrix<double> generatePascalQuad(int order)
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index 3c5dca3e1a84ca1297564b52a7cd312806eda750..0110c5f087efdc4e68740a0e2e531dd7c481d772 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -299,8 +299,6 @@ class polynomialBasis
break;
}
}
- const fullMatrix<double> &getGradientAtFaceIntegrationPoints(int integrationOrder,
- int closureId) const;
static int getTag(int parentTag, int order, bool serendip = false);
};