From cda9793ca10756b8268a70f9c6cea206cc36e1ad Mon Sep 17 00:00:00 2001
From: Thomas De Maet <thomas.demaet@uclouvain.be>
Date: Thu, 8 Mar 2012 10:54:51 +0000
Subject: [PATCH] richards equation stuff + add possibility to compute means
directly on barycenter (for lin elmts)
---
Geo/MElement.cpp | 17 +++++++++++++++++
Geo/MElement.h | 1 +
Geo/MHexahedron.h | 4 ++++
Geo/MLine.h | 4 ++++
Geo/MPoint.h | 4 ++++
Geo/MPrism.h | 4 ++++
Geo/MPyramid.h | 4 ++++
Geo/MQuadrangle.h | 16 ++++++++++++++++
Geo/MTetrahedron.h | 4 ++++
Geo/MTriangle.h | 4 ++++
10 files changed, 62 insertions(+)
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 5881e05096..4c6dbae508 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -158,6 +158,23 @@ SPoint3 MElement::barycenter()
return p;
}
+SPoint3 MElement::barycenterUVW()
+{
+ SPoint3 p(0., 0., 0.);
+ int n = getNumVertices();
+ for(int i = 0; i < n; i++) {
+ double x, y, z;
+ getNode(i, x, y, z);
+ p[0] += x;
+ p[1] += y;
+ p[2] += z;
+ }
+ p[0] /= (double)n;
+ p[1] /= (double)n;
+ p[2] /= (double)n;
+ return p;
+}
+
double MElement::getVolume()
{
int npts;
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 146d116803..4d3cc54a50 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -194,6 +194,7 @@ class MElement
// compute the barycenter
virtual SPoint3 barycenter();
+ virtual SPoint3 barycenterUVW();
// revert the orientation of the element
virtual void revert(){}
diff --git a/Geo/MHexahedron.h b/Geo/MHexahedron.h
index a60ae836a5..9c552acbab 100644
--- a/Geo/MHexahedron.h
+++ b/Geo/MHexahedron.h
@@ -137,6 +137,10 @@ class MHexahedron : public MElement {
default: u = 0.; v = 0.; w = 0.; break;
}
}
+ virtual SPoint3 barycenterUVW()
+ {
+ return SPoint3(0., 0., 0.);
+ }
virtual bool isInside(double u, double v, double w)
{
double tol = _isInsideTolerance;
diff --git a/Geo/MLine.h b/Geo/MLine.h
index 53e77dd72a..a26fbaf268 100644
--- a/Geo/MLine.h
+++ b/Geo/MLine.h
@@ -88,6 +88,10 @@ class MLine : public MElement {
default: u = 0.; break;
}
}
+ virtual SPoint3 barycenterUVW()
+ {
+ return SPoint3(0, 0, 0);
+ }
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
};
diff --git a/Geo/MPoint.h b/Geo/MPoint.h
index e294bd5813..1b09436051 100644
--- a/Geo/MPoint.h
+++ b/Geo/MPoint.h
@@ -46,6 +46,10 @@ class MPoint : public MElement {
{
u = v = w = 0.;
}
+ virtual SPoint3 barycenterUVW()
+ {
+ return SPoint3(0., 0., 0.);
+ }
virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
{
s[0] = 1.;
diff --git a/Geo/MPrism.h b/Geo/MPrism.h
index 33c81597fc..d6e14731ca 100644
--- a/Geo/MPrism.h
+++ b/Geo/MPrism.h
@@ -142,6 +142,10 @@ class MPrism : public MElement {
default: u = 0.; v = 0.; w = 0.; break;
}
}
+ virtual SPoint3 barycenterUVW()
+ {
+ return SPoint3(1/3., 1/3., 0.);
+ }
virtual bool isInside(double u, double v, double w)
{
double tol = _isInsideTolerance;
diff --git a/Geo/MPyramid.h b/Geo/MPyramid.h
index 90676b6f21..065d1186dc 100644
--- a/Geo/MPyramid.h
+++ b/Geo/MPyramid.h
@@ -176,6 +176,10 @@ class MPyramid : public MElement {
default: u = 0.; v = 0.; w = 0.; break;
}
}
+ virtual SPoint3 barycenterUVW()
+ {
+ return SPoint3(0., 0., .2);
+ }
virtual bool isInside(double u, double v, double w)
{
double tol = _isInsideTolerance;
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index a0c934faf3..f3affef3cc 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -128,6 +128,10 @@ void projectInMeanPlane(double *xn, double *yn);
default: u = 0.; v = 0.; break;
}
}
+ virtual SPoint3 barycenterUVW()
+ {
+ return SPoint3(0., 0., 0.);
+ }
virtual void revert()
{
MVertex *tmp = _v[1]; _v[1] = _v[3]; _v[3] = tmp;
@@ -240,6 +244,10 @@ class MQuadrangle8 : public MQuadrangle {
{
num < 4 ? MQuadrangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
}
+ virtual SPoint3 barycenterUVW()
+ {
+ return SPoint3(0., 0., 0.);
+ }
};
/*
@@ -316,6 +324,10 @@ class MQuadrangle9 : public MQuadrangle {
{
num < 4 ? MQuadrangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
}
+ virtual SPoint3 barycenterUVW()
+ {
+ return SPoint3(0., 0., 0.);
+ }
};
/*
@@ -412,6 +424,10 @@ class MQuadrangleN : public MQuadrangle {
{
num < 4 ? MQuadrangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
}
+ virtual SPoint3 barycenterUVW()
+ {
+ return SPoint3(0., 0., 0.);
+ }
};
template <class T>
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index ffae493d34..3f017d8e8f 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -140,6 +140,10 @@ class MTetrahedron : public MElement {
default: u = 0.; v = 0.; w = 0.; break;
}
}
+ virtual SPoint3 barycenterUVW()
+ {
+ return SPoint3(.25, .25, .25);
+ }
virtual bool isInside(double u, double v, double w)
{
double tol = _isInsideTolerance;
diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h
index 2f78ab4fc3..dbcb47e433 100644
--- a/Geo/MTriangle.h
+++ b/Geo/MTriangle.h
@@ -136,6 +136,10 @@ class MTriangle : public MElement {
default: u = 0.; v = 0.; break;
}
}
+ virtual SPoint3 barycenterUVW()
+ {
+ return SPoint3(1/3., 1/3., 0.);
+ }
virtual bool isInside(double u, double v, double w)
{
double tol = _isInsideTolerance;
--
GitLab