Added ComputeMetricTensor and GetMetricTensor for the g_ij (as opposed

to GetGmat for the g^ij). Added creation of the manifold mass and
derivative matrices.

library/LocalRegions/Expansion2D.cpp

@@ -1091,23 +1091,66 @@ namespace Nektar
 			//mass matrix M[i,j] = (phi_i, g_11 phi_j)
 			case StdRegions::eManifoldMass00:
 				{
-					Array<TwoD, const NekDouble> gmat = m_metricinfo->GetGmat(GetPointsKeys());//!!this is inverted metric tensor, need to write a function to compute g_ij one
-					StdRegions::VarCoeffsMap vcMap;
-					vcMap[StdRegions::eVarCoeffMass] = gmat[0][0];
+					int nqtot = GetTotPoints();
+					Array<TwoD, const NekDouble> gmat = m_metricinfo->GetMetricTensor(GetPointsKeys());
+					//passing metric info through the VarCoeffMap in matrix key
+					StdRegions::VarCoeffMap vcMap;
+					vcMap[StdRegions::eVarCoeffMetric] = Array<OneD, NekDouble>(nqtot); 
+					Vmath::Vcopy(nqtot, &gmat[0][0], 1, &vcMap[StdRegions::eVarCoeffMetric][0], 1);
 					StdRegions::StdMatrixKey masskey(StdRegions::eMass, DetShapeType(), *this, StdRegions::NullConstFactorMap, vcMap);
-					returnval = StdRegions::CreateGeneralMatrix(masskey);
+					returnval = StdExpansion::CreateGeneralMatrix(masskey);
 				}
 				break;
 			//mass matrix M[i,j] = (phi_i, g_21 phi_j). Keep in mind that g_12 = g_21 due to the symmetry of the metric tensor
 			case StdRegions::eManifoldMass10:
+				{
+					int nqtot = GetTotPoints();
+					Array<TwoD, const NekDouble> gmat = m_metricinfo->GetMetricTensor(GetPointsKeys());
+					//passing metric info through the VarCoeffMap in matrix key
+					StdRegions::VarCoeffMap vcMap;
+					vcMap[StdRegions::eVarCoeffMetric] = Array<OneD, NekDouble>(nqtot); 
+					Vmath::Vcopy(nqtot, &gmat[1][0], 1, &vcMap[StdRegions::eVarCoeffMetric][0], 1);
+					StdRegions::StdMatrixKey masskey(StdRegions::eMass, DetShapeType(), *this, StdRegions::NullConstFactorMap, vcMap);
+					returnval = StdExpansion::CreateGeneralMatrix(masskey);
+				}
 				break;
 			//mass matrix M[i,j] = (phi_i, g_22 phi_j)
 			case StdRegions::eManifoldMass11:
+				{
+					int nqtot = GetTotPoints();
+					Array<TwoD, const NekDouble> gmat = m_metricinfo->GetMetricTensor(GetPointsKeys());
+					//passing metric info through the VarCoeffMap in matrix key
+					StdRegions::VarCoeffMap vcMap;
+					vcMap[StdRegions::eVarCoeffMetric] = Array<OneD, NekDouble>(nqtot); 
+					Vmath::Vcopy(nqtot, &gmat[4][0], 1, &vcMap[StdRegions::eVarCoeffMetric][0], 1);
+					StdRegions::StdMatrixKey masskey(StdRegions::eMass, DetShapeType(), *this, StdRegions::NullConstFactorMap, vcMap);
+					returnval = StdExpansion::CreateGeneralMatrix(masskey);
+				}
 				break;
 			//weak derivative matrices for manifod D_k[i,j] = (phi_i, d( sqrt(g) phi_j )/d xi_k)
 			case StdRegions::eManifoldWeakDeriv0:
+				{
+					int nqtot = GetTotPoints();
+					Array<OneD, const NekDouble> jac = m_metricinfo->GetJac(GetPointsKeys());
+					//passing metric info through the VarCoeffMap in matrix key
+					StdRegions::VarCoeffMap vcMap;
+					vcMap[StdRegions::eVarCoeffMetric] = Array<OneD, NekDouble>(nqtot); 
+					Vmath::Vcopy(nqtot, &jac[0], 1, &vcMap[StdRegions::eVarCoeffMetric][0], 1);
+					StdRegions::StdMatrixKey weakDkey(StdRegions::eWeakDeriv0, DetShapeType(), *this, StdRegions::NullConstFactorMap, vcMap);
+					returnval = StdExpansion::CreateGeneralMatrix(weakDkey);
+				}
 				break;
 			case StdRegions::eManifoldWeakDeriv1:
+				{
+					int nqtot = GetTotPoints();
+					Array<OneD, const NekDouble> jac = m_metricinfo->GetJac(GetPointsKeys());
+					//passing metric info through the VarCoeffMap in matrix key
+					StdRegions::VarCoeffMap vcMap;
+					vcMap[StdRegions::eVarCoeffMetric] = Array<OneD, NekDouble>(nqtot); 
+					Vmath::Vcopy(nqtot, &jac[0], 1, &vcMap[StdRegions::eVarCoeffMetric][0], 1);
+					StdRegions::StdMatrixKey weakDkey(StdRegions::eWeakDeriv0, DetShapeType(), *this, StdRegions::NullConstFactorMap, vcMap);
+					returnval = StdExpansion::CreateGeneralMatrix(weakDkey);
+				}
 				break;
             default:
                 ASSERTL0(false,"This matrix type cannot be generated from this class");

library/SpatialDomains/GeomFactors.cpp

@@ -307,6 +307,71 @@ namespace Nektar
             return jac;
         }
 
+        /**
+         * This routine returns a two-dimensional array of values specifying
+         * the metric terms associated with the coordinate mapping of
+         * the corresponding reference region to the physical element. These
+         * terms correspond to the \f$g_{ij}\f$ terms in \cite CaYaKiPeSh13 and,
+         * in the case of an embedded manifold, map contravariant quantities to
+         * covariant quantities. The leading index of the array is the index
+         * of the term in the tensor numbered as
+         * \f[\left(\begin{array}{ccc}
+         *    0 & 1 & 2 \\
+         *    1 & 3 & 4 \\
+         *    2 & 4 & 5
+         * \end{array}\right)\f].
+         * The second dimension is either of size 1 in the case of elements
+         * having #GeomType #eRegular, or of size equal to the number of
+         * quadrature points for #eDeformed elements.
+         *
+         * @see [Wikipedia "Covariance and Contravariance of Vectors"]
+         *      (http://en.wikipedia.org/wiki/Covariance_and_contravariance_of_vectors)
+         * @returns             Two-dimensional array containing the inverse
+         *                      metric tensor of the coordinate mapping.
+         */
+        Array<TwoD, NekDouble> GeomFactors::ComputeMetricTensor(
+                const LibUtilities::PointsKeyVector &keyTgt) const
+        {
+            ASSERTL1(keyTgt.size() == m_expDim,
+                     "Dimension of target point distribution does not match "
+                     "expansion dimension.");
+
+            int i = 0, j = 0, k = 0, l = 0;
+            int ptsTgt   = 1;
+
+            if (m_type == eDeformed)
+            {
+                // Allocate storage and compute number of points
+                for (i = 0; i < m_expDim; ++i)
+                {
+                    ptsTgt   *= keyTgt[i].GetNumPoints();
+                }
+            }
+
+            // Get derivative at geometry points
+            DerivStorage deriv = ComputeDeriv(keyTgt);
+
+            Array<TwoD, NekDouble> gmat(m_expDim*m_expDim, ptsTgt, 0.0);
+
+            // Compute g_{ij} as t_i \cdot t_j and store in tmp
+            for (i = 0, l = 0; i < m_expDim; ++i)
+            {
+                for (j = 0; j < m_expDim; ++j, ++l)
+                {
+                    for (k = 0; k < m_coordDim; ++k)
+                    {
+                        Vmath::Vvtvp(ptsTgt, &deriv[i][k][0], 1,
+                                             &deriv[j][k][0], 1,
+                                             &gmat[l][0],      1,
+                                             &gmat[l][0],      1);
+                    }
+                }
+            }
+
+            return gmat;
+        }
+
+
 
         /**
          * This routine returns a two-dimensional array of values specifying

library/SpatialDomains/GeomFactors.h

@@ -107,7 +107,11 @@ namespace SpatialDomains
             inline const Array<OneD, const NekDouble> GetJac(
                     const LibUtilities::PointsKeyVector &keyTgt);
 
-            /// Return the Laplacian coefficients \f$g_{ij}\f$.
+            /// Return the metric tensor \f$g_{ij}\f$.
+            inline const Array<TwoD, const NekDouble> GetMetricTensor(
+                    const LibUtilities::PointsKeyVector &keyTgt);
+
+            /// Return the Laplacian coefficients (inverse metric tensor) \f$g^{ij}\f$.
             inline const Array<TwoD, const NekDouble> GetGmat(
                     const LibUtilities::PointsKeyVector &keyTgt);
 
@@ -159,7 +163,11 @@ namespace SpatialDomains
             SPATIAL_DOMAINS_EXPORT Array<OneD, NekDouble> ComputeJac(
                     const LibUtilities::PointsKeyVector &keyTgt) const;
 
-            /// Computes the Laplacian coefficients \f$g_{ij}\f$.
+            /// Computes the metric tensor \f$g_{ij}\f$.
+            SPATIAL_DOMAINS_EXPORT Array<TwoD, NekDouble> ComputeMetricTensor(
+                    const LibUtilities::PointsKeyVector &keyTgt) const;
+
+            /// Computes the Laplacian coefficients (inverse metric tensor) \f$g^{ij}\f$.
             SPATIAL_DOMAINS_EXPORT Array<TwoD, NekDouble> ComputeGmat(
                     const LibUtilities::PointsKeyVector &keyTgt) const;
 
@@ -235,6 +243,18 @@ namespace SpatialDomains
     }
 
 
+    /**
+     * @param   keyTgt      Target point distributions.
+     * @returns             Metric tensor evaluated at target point
+     *                      distributions.
+     * @see                 GeomFactors::ComputeMetricTensor
+     */
+    inline const Array<TwoD, const NekDouble> GeomFactors::GetMetricTensor(
+            const LibUtilities::PointsKeyVector &keyTgt)
+    {
+        return ComputeMetricTensor(keyTgt);
+    }
+
     /**
      * @param   keyTgt      Target point distributions.
      * @returns             Inverse metric tensor evaluated at target point

library/StdRegions/StdRegions.hpp

@@ -209,7 +209,8 @@ namespace Nektar
             eVarCoeffD02,
             eVarCoeffD12,
             eVarCoeffVelX,
-            eVarCoeffVelY
+            eVarCoeffVelY,
+			eVarCoeffMetric
         };
 
         const char* const VarCoeffTypeMap[] = {
@@ -223,7 +224,8 @@ namespace Nektar
             "VarCoeffD02",
             "VarCoeffD12",
             "VarCoeffVelX",
-            "VarCoeffVelY"
+            "VarCoeffVelY",
+			"VarCoeffMetric"
         };
         typedef std::map<StdRegions::VarCoeffType, Array<OneD, NekDouble> > VarCoeffMap;
         static VarCoeffMap NullVarCoeffMap;