From 71f01c1c40f14a0767e544ae55779898c08c8138 Mon Sep 17 00:00:00 2001
From: Thomas De Maet <thomas.demaet@uclouvain.be>
Date: Fri, 30 Mar 2012 14:49:23 +0000
Subject: [PATCH] Ground-surface water coupling
---
Geo/MLine.cpp | 5 +++++
Geo/MLine.h | 1 +
Geo/MQuadrangle.cpp | 15 ++++++++++++++-
Geo/MQuadrangle.h | 3 ++-
Geo/MTriangle.cpp | 11 +++++++++++
Geo/MTriangle.h | 1 +
6 files changed, 34 insertions(+), 2 deletions(-)
diff --git a/Geo/MLine.cpp b/Geo/MLine.cpp
index 5091bc18cc..3ffe7f1e10 100644
--- a/Geo/MLine.cpp
+++ b/Geo/MLine.cpp
@@ -65,3 +65,8 @@ double MLine::getLength()
{
return _v[0]->distance(_v[1]);
}
+
+double MLine::getVolume()
+{
+ return getLength();
+}
diff --git a/Geo/MLine.h b/Geo/MLine.h
index a26fbaf268..8e60f50a6a 100644
--- a/Geo/MLine.h
+++ b/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;
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index c448fabab4..6f41fe2dc3 100644
--- a/Geo/MQuadrangle.cpp
+++ b/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)
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index f3affef3cc..4484c1dad8 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -38,7 +38,7 @@ class MQuadrangle : public MElement {
v[2] = _v[2];
v[3] = _v[3];
}
-void projectInMeanPlane(double *xn, double *yn);
+ void projectInMeanPlane(double *xn, double *yn);
public :
MQuadrangle(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3, int num=0, int part=0)
: MElement(num, part)
@@ -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;
diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp
index a15d74f8ef..98f977d4f1 100644
--- a/Geo/MTriangle.cpp
+++ b/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)
diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h
index dbcb47e433..54a298904e 100644
--- a/Geo/MTriangle.h
+++ b/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] = {
--
GitLab