Ground-surface water coupling

Geo/MLine.cpp

@@ -65,3 +65,8 @@ double MLine::getLength()
 {
   return _v[0]->distance(_v[1]);
 }
+
+double MLine::getVolume()
+{
+  return getLength();
+}

Geo/MLine.h

@@ -39,6 +39,7 @@ class MLine : public MElement {
   virtual MVertex *getVertex(int num){ return _v[num]; }
   virtual double getInnerRadius(); // half-length of segment line
   virtual double getLength(); // length of segment line
+  virtual double getVolume();
   virtual void getVertexInfo(const MVertex * vertex, int &ithVertex) const
   {
     ithVertex = _v[0] == vertex ? 0 : 1;

Geo/MQuadrangle.cpp

@@ -94,7 +94,20 @@ int MQuadrangleN::getNumEdgesRep(){ return 4 * CTX::instance()->mesh.numSubEdges
 int MQuadrangle8::getNumEdgesRep(){ return 4 * CTX::instance()->mesh.numSubEdges; }
 int MQuadrangle9::getNumEdgesRep(){ return 4 * CTX::instance()->mesh.numSubEdges; }
 
-
+double MQuadrangle::getVolume()
+{
+  if(getNumVertices() > 4)
+    return MElement::getVolume();
+  double a = _v[0]->distance(_v[1]);
+  double b = _v[1]->distance(_v[2]);
+  double c = _v[2]->distance(_v[3]);
+  double d = _v[3]->distance(_v[0]);
+  double m = _v[0]->distance(_v[2]);
+  double n = _v[1]->distance(_v[3]);
+  double mn = 2. * m * n;
+  double abcd = a*a - b*b + c*c - d*d;
+  return sqrt( mn*mn - abcd*abcd ) / 4.;
+}
 
 static void _myGetEdgeRep(MQuadrangle *q, int num, double *x, double *y, double *z,
                           SVector3 *n, int numSubEdges)

Geo/MQuadrangle.h

@@ -132,6 +132,7 @@ void projectInMeanPlane(double *xn, double *yn);
   {
     return SPoint3(0., 0., 0.);
   }
+  virtual double getVolume();
   virtual void revert()
   {
     MVertex *tmp = _v[1]; _v[1] = _v[3]; _v[3] = tmp;

Geo/MTriangle.cpp

@@ -25,6 +25,17 @@ SPoint3 MTriangle::circumcenter()
   return SPoint3(res[0], res[1], res[2]);
 }
 
+double MTriangle::getVolume()
+{
+  if(getNumVertices() > 3)
+    return MElement::getVolume();
+  SPoint3 p0(_v[0]->x(), _v[0]->y(), _v[0]->z());
+  SPoint3 p1(_v[1]->x(), _v[1]->y(), _v[1]->z());
+  SPoint3 p2(_v[2]->x(), _v[2]->y(), _v[2]->z());
+  SVector3 v1(p0, p1), v2(p0, p2);
+  return norm(crossprod(v1, v2)) / 2.;
+}
+
 double MTriangle::distoShapeMeasure()
 {
 #if defined(HAVE_MESH)

Geo/MTriangle.h

@@ -149,6 +149,7 @@ class MTriangle : public MElement {
   }
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
   virtual SPoint3 circumcenter();
+  virtual double getVolume();
   static int edges_tri(const int edge, const int vert)
   {
     static const int e[3][2] = {