dg : remove all constant matrices from dgGroupOf...

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];

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]);

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)

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);
 };