MElement.cpp

  case MSH_QUA_8:  return new MQuadrangle8(v, num, part);
  case MSH_QUA_9:  return new MQuadrangle9(v, num, part);
  case MSH_QUA_12: return new MQuadrangleN(v, 3, num, part);
  case MSH_QUA_16: return new MQuadrangleN(v, 3, num, part);
  case MSH_QUA_25: return new MQuadrangleN(v, 4, num, part);
  case MSH_QUA_36: return new MQuadrangleN(v, 5, num, part);
  case MSH_QUA_49: return new MQuadrangleN(v, 6, num, part);
  case MSH_QUA_64: return new MQuadrangleN(v, 7, num, part);
  case MSH_QUA_81: return new MQuadrangleN(v, 8, num, part);
  case MSH_QUA_100:return new MQuadrangleN(v, 9, num, part);
  case MSH_QUA_121:return new MQuadrangleN(v, 10, num, part);
  case MSH_POLYG_: return new MPolygon(v, num, part, owner, parent);
  case MSH_POLYG_B:return new MPolygonBorder(v, num, part, d1, d2);
  case MSH_TET_4:  return new MTetrahedron(v, num, part);
  case MSH_TET_10: return new MTetrahedron10(v, num, part);
  case MSH_HEX_8:  return new MHexahedron(v, num, part);
  case MSH_HEX_20: return new MHexahedron20(v, num, part);
  case MSH_HEX_27: return new MHexahedron27(v, num, part);
  case MSH_PRI_6:  return new MPrism(v, num, part);
  case MSH_PRI_15: return new MPrism15(v, num, part);
  case MSH_PRI_18: return new MPrism18(v, num, part);
  case MSH_PYR_5:  return new MPyramid(v, num, part);
  case MSH_PYR_13: return new MPyramid13(v, num, part);
  case MSH_PYR_14: return new MPyramid14(v, num, part);
  case MSH_TET_20: return new MTetrahedronN(v, 3, num, part);
  case MSH_TET_34: return new MTetrahedronN(v, 3, num, part);
  case MSH_TET_35: return new MTetrahedronN(v, 4, num, part);
  case MSH_TET_52: return new MTetrahedronN(v, 5, num, part);
  case MSH_TET_56: return new MTetrahedronN(v, 5, num, part);
  case MSH_TET_84: return new MTetrahedronN(v, 6, num, part);
  case MSH_TET_120: return new MTetrahedronN(v, 7, num, part);
  case MSH_TET_165: return new MTetrahedronN(v, 8, num, part);
  case MSH_TET_220: return new MTetrahedronN(v, 9, num, part);
  case MSH_TET_286: return new MTetrahedronN(v, 10, num, part);
  case MSH_POLYH_: return new MPolyhedron(v, num, part, owner, parent);
  default:         return 0;
  }
}

/*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);
  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 (i, j) = f[j];
    }
  }
  return mat;
}*/

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

#include "Bindings.h"

void MElement::registerBindings(binding *b)
{
  classBinding *cb = b->addClass<MElement>("MElement");
  cb->setDescription("A mesh element.");
  methodBinding *cm;
  cm = cb->addMethod("getNum",&MElement::getNum);
  cm->setDescription("return the tag of the element");
  cm = cb->addMethod("getNumVertices", &MElement::getNumVertices);
  cm->setDescription("get the number of vertices of this element");
  cm = cb->addMethod("getVertex", &MElement::getVertex);
  cm->setDescription("return the i-th vertex of this element");
  cm->setArgNames("i",NULL);
  cm = cb->addMethod("getType", &MElement::getType);
  cm->setDescription("get the type of the element");
  cm = cb->addMethod("getTypeForMSH", &MElement::getTypeForMSH);
  cm->setDescription("get the gmsh type of the element");
  cm = cb->addMethod("getPartition", &MElement::getPartition);
  cm->setDescription("get the partition to which the element belongs");
  cm = cb->addMethod("setPartition", &MElement::setPartition);
  cm->setDescription("set the partition to which the element belongs");
  cm->setArgNames("iPartition",NULL);
  cm = cb->addMethod("getPolynomialOrder", &MElement::getPolynomialOrder);
  cm->setDescription("return the polynomial order the element");
  cm = cb->addMethod("getDim", &MElement::getDim);
  cm->setDescription("return the geometrical dimension of the element");
  cm = cb->addMethod("getJacobianDeterminant", &MElement::getJacobianDeterminant);
  cm->setDescription("return the jacobian of the determinant of the transformation");
  cm->setArgNames("u","v","w",NULL);
}