updated msubelement

Geo/MSubElement.cpp

@@ -12,8 +12,7 @@
 
 MSubTetrahedron::~MSubTetrahedron()
 {
-//   if(_owner)
-//     delete _orig;
+  if(_intpt) delete [] _intpt;
 }
 
 const nodalBasis* MSubTetrahedron::getFunctionSpace(int order) const
@@ -50,17 +49,20 @@ void MSubTetrahedron::getThirdDerivativeShapeFunctions(double u, double v, doubl
 
 int MSubTetrahedron::getNumShapeFunctions()
 {
-  if(_orig) _orig->getNumShapeFunctions();
+  if(_orig) return _orig->getNumShapeFunctions();
+  return 0;
 }
 
 int MSubTetrahedron::getNumPrimaryShapeFunctions()
 {
-  if(_orig) _orig->getNumPrimaryShapeFunctions();
+  if(_orig) return _orig->getNumPrimaryShapeFunctions();
+  return 0;
 }
 
 MVertex* MSubTetrahedron::getShapeFunctionNode(int i)
 {
   if(_orig) return _orig->getShapeFunctionNode(i);
+  return 0;
 }
 
 bool MSubTetrahedron::isInside(double u, double v, double w)
@@ -95,6 +97,8 @@ void MSubTetrahedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
   int nptsi;
   IntPt *ptsi;
 
+  // work in the parametric space of the parent element
+  // (i) get the parametric coordinates of MSubElement vertices in this space
   double v_uvw[4][3];
   for(int i=0; i<4; ++i)
   {
@@ -102,19 +106,20 @@ void MSubTetrahedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
     double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
     _orig->xyz2uvw(v_xyz, v_uvw[i]);
   }
+  // (ii) build a MElement on the vertices with these coordinates in this space
   MVertex v0(v_uvw[0][0], v_uvw[0][1], v_uvw[0][2]);
   MVertex v1(v_uvw[1][0], v_uvw[1][1], v_uvw[1][2]);
   MVertex v2(v_uvw[2][0], v_uvw[2][1], v_uvw[2][2]);
   MVertex v3(v_uvw[3][0], v_uvw[3][1], v_uvw[3][2]);
   MTetrahedron tt(&v0, &v1, &v2, &v3);
+  // (iii) get the integration points of the new element in its parametric space 
   tt.getIntegrationPoints(pOrder, &nptsi, &ptsi);
-  double jac[3][3];
+  // (iv) get the coordinates of these points in the parametric space of parent element
   for(int ip=0; ip<nptsi; ++ip)
   {
     const double u = ptsi[ip].pt[0];
     const double v = ptsi[ip].pt[1];
     const double w = ptsi[ip].pt[2];
-    tt.getJacobian(u, v, w, jac);
     SPoint3 p; tt.pnt(u, v, w, p);
     _intpt[*npts + ip].pt[0] = p.x();
     _intpt[*npts + ip].pt[1] = p.y();
@@ -130,8 +135,7 @@ void MSubTetrahedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
 
 MSubTriangle::~MSubTriangle()
 {
-//   if(_owner)
-//     delete _orig;
+  if(_intpt) delete [] _intpt;
 }
 
 const nodalBasis* MSubTriangle::getFunctionSpace(int order) const
@@ -168,17 +172,20 @@ void MSubTriangle::getThirdDerivativeShapeFunctions(double u, double v, double w
 
 int MSubTriangle::getNumShapeFunctions()
 {
-  if(_orig) _orig->getNumShapeFunctions();
+  if(_orig) return _orig->getNumShapeFunctions();
+  return 0;
 }
 
 int MSubTriangle::getNumPrimaryShapeFunctions()
 {
-  if(_orig) _orig->getNumPrimaryShapeFunctions();
+  if(_orig) return _orig->getNumPrimaryShapeFunctions();
+  return 0;
 }
 
 MVertex* MSubTriangle::getShapeFunctionNode(int i)
 {
   if(_orig) return _orig->getShapeFunctionNode(i);
+  return 0;
 }
 
 bool MSubTriangle::isInside(double u, double v, double w)
@@ -212,6 +219,8 @@ void MSubTriangle::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
   int nptsi;
   IntPt *ptsi;
 
+  // work in the parametric space of the parent element
+  // (i) get the parametric coordinates of MSubElement vertices in this space
   double v_uvw[3][3];
   for(int i=0; i<3; ++i)
   {
@@ -219,18 +228,19 @@ void MSubTriangle::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
     double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
     _orig->xyz2uvw(v_xyz, v_uvw[i]);
   }
+  // (ii) build a MElement on the vertices with these coordinates in this space
   MVertex v0(v_uvw[0][0], v_uvw[0][1], v_uvw[0][2]);
   MVertex v1(v_uvw[1][0], v_uvw[1][1], v_uvw[1][2]);
   MVertex v2(v_uvw[2][0], v_uvw[2][1], v_uvw[2][2]);
   MTriangle t(&v0, &v1, &v2);
+  // (iii) get the integration points of the new element in its parametric space 
   t.getIntegrationPoints(pOrder, &nptsi, &ptsi);
-  double jac[3][3];
+  // (iv) get the coordinates of these points in the parametric space of parent element
   for(int ip=0; ip<nptsi; ++ip)
   {
     const double u = ptsi[ip].pt[0];
     const double v = ptsi[ip].pt[1];
     const double w = ptsi[ip].pt[2];
-    t.getJacobian(u, v, w, jac);
     SPoint3 p; t.pnt(u, v, w, p);
     _intpt[*npts + ip].pt[0] = p.x();
     _intpt[*npts + ip].pt[1] = p.y();
@@ -245,8 +255,7 @@ void MSubTriangle::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
 
 MSubLine::~MSubLine()
 {
-//   if(_owner)
-//     delete _orig;
+  if(_intpt) delete [] _intpt;
 }
 
 const nodalBasis* MSubLine::getFunctionSpace(int order) const
@@ -283,17 +292,20 @@ void MSubLine::getThirdDerivativeShapeFunctions(double u, double v, double w, do
 
 int MSubLine::getNumShapeFunctions()
 {
-  if(_orig) _orig->getNumShapeFunctions();
+  if(_orig) return _orig->getNumShapeFunctions();
+  return 0;
 }
 
 int MSubLine::getNumPrimaryShapeFunctions()
 {
-  if(_orig) _orig->getNumPrimaryShapeFunctions();
+  if(_orig) return _orig->getNumPrimaryShapeFunctions();
+  return 0;
 }
 
 MVertex* MSubLine::getShapeFunctionNode(int i)
 {
   if(_orig) return _orig->getShapeFunctionNode(i);
+  return 0;
 }
 
 bool MSubLine::isInside(double u, double v, double w)
@@ -325,6 +337,9 @@ void MSubLine::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
   _intpt = new IntPt[getNGQLPts(pOrder)];
   int nptsi;
   IntPt *ptsi;
+
+  // work in the parametric space of the parent element
+  // (i) get the parametric coordinates of MSubElement vertices in this space
   double v_uvw[2][3];
   for(int i=0; i<2; ++i)
   {
@@ -332,10 +347,13 @@ void MSubLine::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
     double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
     _orig->xyz2uvw(v_xyz, v_uvw[i]);
   }
+  // (ii) build a MElement on the vertices with these coordinates in this space
   MVertex v0(v_uvw[0][0], v_uvw[0][1], v_uvw[0][2]);
   MVertex v1(v_uvw[1][0], v_uvw[1][1], v_uvw[1][2]);
   MLine l(&v0, &v1);
+  // (iii) get the integration points of the new element in its parametric space 
   l.getIntegrationPoints(pOrder, &nptsi, &ptsi);
+  // (iv) get the coordinates of these points in the parametric space of parent element
   for(int ip=0; ip<nptsi; ++ip)
   {
     const double u = ptsi[ip].pt[0];
@@ -356,8 +374,7 @@ void MSubLine::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
 
 MSubPoint::~MSubPoint()
 {
-//   if(_owner)
-//     delete _orig;
+  if(_intpt) delete [] _intpt;
 }
 
 const nodalBasis* MSubPoint::getFunctionSpace(int order) const
@@ -394,17 +411,20 @@ void MSubPoint::getThirdDerivativeShapeFunctions(double u, double v, double w, d
 
 int MSubPoint::getNumShapeFunctions()
 {
-  if(_orig) _orig->getNumShapeFunctions();
+  if(_orig) return _orig->getNumShapeFunctions();
+  return 0;
 }
 
 int MSubPoint::getNumPrimaryShapeFunctions()
 {
-  if(_orig) _orig->getNumPrimaryShapeFunctions();
+  if(_orig) return _orig->getNumPrimaryShapeFunctions();
+  return 0;
 }
 
 MVertex* MSubPoint::getShapeFunctionNode(int i)
 {
   if(_orig) return _orig->getShapeFunctionNode(i);
+  return 0;
 }
 
 bool MSubPoint::isInside(double u, double v, double w)
@@ -425,11 +445,21 @@ bool MSubPoint::isInside(double u, double v, double w)
 
 void MSubPoint::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
 {
-  static IntPt GQL[1];
-  GQL[0].pt[0] = 0;
-  GQL[0].pt[1] = 0;
-  GQL[0].pt[2] = 0;
-  GQL[0].weight = 1;
+  // invariable regardless of the order
+  if(!_intpt)
+  {
+    if(!_orig) return;
+    _intpt = new IntPt[1];
+    // work in the parametric space of the parent element
+    MVertex *vi = getVertex(0);
+    double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
+    double v_uvw[3];
+    _orig->xyz2uvw(v_xyz, v_uvw);
+    _intpt[0].pt[0] = v_uvw[0];
+    _intpt[0].pt[1] = v_uvw[1];
+    _intpt[0].pt[2] = v_uvw[2];
+    _intpt[0].weight = 1;
+  }
   *npts = 1;
-  *pts = GQL;
+  *pts = _intpt;
 }