added raviart-thomas for interval

src/Discretization.C

@@ -311,6 +311,8 @@ void Discretization::fill(int m, int dim) {
         S[1][1] = Get_Shape("RT0quad", "Qquad4", "Qint1");
         S[4][0] = Get_Shape("P0Hex", "Qhex8", "Qint1");
         S[4][1] = Get_Shape("RT0hex", "Qhex8", "Qint1");
+        S[5][0] = Get_Shape("P0Int", "Qint2", "Qpoint1");
+        S[5][1] = Get_Shape("RT0int", "Qint2", "Qpoint1");
     } else if (name == "VectorP1RT0") {
         S[0][0] = Get_Shape("P1Tri", "Qtri3", "Qint1");
         S[0][1] = Get_Shape("RT0tri", "Qtri1", "Qint1");

src/shapes/RaviartThomasShape.cpp

+#include "RaviartThomasShape.hpp"
+
+double RT0tri::LocalDiv(const Point &z, int i) const {
+    return 2.0;
+}
+
+VectorField RT0quad::LocalVector(const Point &z, int i) const {
+    switch (i) {
+        case 0:return VectorField(0.0, z[1] - 1.0);
+        case 1:return VectorField(z[0], 0.0);
+        case 2:return VectorField(0.0, z[1]);
+        case 3:return VectorField(z[0] - 1.0, 0.0);
+    }
+    Exit("Not implemented");
+    return 0.0;
+}
+
+double RT0quad::LocalDiv(const Point &z, int i) const {
+    return 1.0;
+}
+
+VectorField RT0tet::LocalVector(const Point &z, int i) const {
+    return zero;
+}
+
+double RT0tet::LocalDiv(const Point &z, int i) const {
+    return 0.0;
+}
+
+VectorField RT0hex::LocalVector(const Point &z, int i) const {
+    switch (i) {
+        case 0:return VectorField(0.0, 0.0, z[2] - 1.0);
+        case 1:return VectorField(0.0, z[1] - 1.0, 0.0);
+        case 2:return VectorField(z[0], 0.0, 0.0);
+        case 3:return VectorField(0.0, z[1], 0.0);
+        case 4:return VectorField(z[0] - 1.0, 0.0, 0.0);
+        case 5:return VectorField(0.0, 0.0, z[2]);
+    }
+    Exit("Not implemented");
+    return 0.0;
+}
+
+VectorField RT0tri::LocalVector(const Point &z, int i) const {
+    switch (i) {
+        case 0:return VectorField(z[0], z[1] - 1.0);
+        case 1:return VectorField(z[0], z[1]);
+        case 2:return VectorField(z[0] - 1.0, z[1]);
+    }
+    Exit("Not implemented");
+    return 0.0;
+}
\ No newline at end of file

src/shapes/RaviartThomasShape.hpp

+#ifndef RAVIARTTHOMASSHAPE_HPP
+#define RAVIARTTHOMASSHAPE_HPP
+
+#include "Shape.h"
+
+
+class RT0Int : public Shape {
+public:
+    RT0Int(const Quadrature &quadrature, const Quadrature &faceQuadrature) :
+        Shape(quadrature, faceQuadrature, 1, 0, 0) {
+    }
+};
+
+class RT0tri : public Shape {
+public:
+    const char *Name() const { return "RT0tri"; }
+
+    RT0tri(const Quadrature &Q, const Quadrature &FaceQ = GetQuadrature("Qint1"))
+        :
+        Shape(Q, FaceQ, 2, 3, 3) {
+        // fill(*this);
+    }
+
+    VectorField LocalVector(const Point &z, int i) const;
+
+    double LocalDiv(const Point &z, int i) const;
+};
+
+class RT0quad : public Shape {
+public:
+    const char *Name() const { return "RT0quad"; }
+
+    RT0quad(const Quadrature &Q,
+            const Quadrature &FaceQ = GetQuadrature("Qint1"))
+        :
+        Shape(Q, FaceQ, 2, 4, 4) {
+        // fill(*this);
+    }
+
+    VectorField LocalVector(const Point &z, int i) const;
+
+    double LocalDiv(const Point &z, int i) const;
+};
+
+class RT0tet : public Shape {
+public:
+    const char *Name() const { return "RT0tet"; }
+
+    RT0tet(const Quadrature &Q,
+           const Quadrature &FaceQ = GetQuadrature("Qint1"))
+        :
+        Shape(Q, FaceQ, 3, 4, 4) {
+        // fill(*this);
+    }
+
+    VectorField LocalVector(const Point &z, int i) const;
+
+    double LocalDiv(const Point &z, int i) const;
+};
+
+class RT0hex : public Shape {
+public:
+    const char *Name() const { return "RT0hex"; }
+
+    RT0hex(const Quadrature &Q,
+           const Quadrature &FaceQ = GetQuadrature("Qint1"))
+        :
+        Shape(Q, FaceQ, 3, 6, 6) {
+        // fill(*this);
+    }
+
+    VectorField LocalVector(const Point &z, int i) const;
+
+    double LocalDiv(const Point &z, int i) const;
+};
+
+#endif //RAVIARTTHOMASSHAPE_HPP

src/shapes/Shape.C

@@ -7,6 +7,7 @@
 #include "HexahedronShape.h"
 #include "TriangularShape.h"
 #include "TetrahedronShape.h"
+#include "RaviartThomasShape.hpp"
 
 
 Point Shape::Eval(const Point &z, const vector<Point> &x) const {
@@ -825,70 +826,6 @@ public:
     }
 };
 
-class RT0tri: public Shape {
-public:
-    const char *Name() const { return "RT0tri"; }
-
-    RT0tri(const Quadrature &Q,
-           const Quadrature &FaceQ = GetQuadrature("Qint1"))
-        :
-        Shape(Q, FaceQ, 2, 3, 3) {
-        // fill(*this);
-    }
-
-    VectorField LocalVector(const Point &z, int i) const;
-
-    double LocalDiv(const Point &z, int i) const;
-};
-
-class RT0quad: public Shape {
-public:
-    const char *Name() const { return "RT0quad"; }
-
-    RT0quad(const Quadrature &Q,
-            const Quadrature &FaceQ = GetQuadrature("Qint1"))
-        :
-        Shape(Q, FaceQ, 2, 4, 4) {
-        // fill(*this);
-    }
-
-    VectorField LocalVector(const Point &z, int i) const;
-
-    double LocalDiv(const Point &z, int i) const;
-};
-
-class RT0tet: public Shape {
-public:
-    const char *Name() const { return "RT0tet"; }
-
-    RT0tet(const Quadrature &Q,
-           const Quadrature &FaceQ = GetQuadrature("Qint1"))
-        :
-        Shape(Q, FaceQ, 3, 4, 4) {
-        // fill(*this);
-    }
-
-    VectorField LocalVector(const Point &z, int i) const;
-
-    double LocalDiv(const Point &z, int i) const;
-};
-
-class RT0hex: public Shape {
-public:
-    const char *Name() const { return "RT0hex"; }
-
-    RT0hex(const Quadrature &Q,
-           const Quadrature &FaceQ = GetQuadrature("Qint1"))
-        :
-        Shape(Q, FaceQ, 3, 6, 6) {
-        // fill(*this);
-    }
-
-    VectorField LocalVector(const Point &z, int i) const;
-
-    double LocalDiv(const Point &z, int i) const;
-};
-
 VectorField P2curl_sp_hex::LocalGradient(const Point &z, int i) const {
     switch (i) {
         case 0:
@@ -1477,56 +1414,6 @@ VectorField P2curl_tet::LocalVector(const Point &z, int i) const {
     return zero;
 }
 
-VectorField RT0tri::LocalVector(const Point &z, int i) const {
-    switch (i) {
-        case 0:return VectorField(z[0], z[1] - 1.0);
-        case 1:return VectorField(z[0], z[1]);
-        case 2:return VectorField(z[0] - 1.0, z[1]);
-    }
-    Exit("Not implemented");
-    return 0.0;
-}
-
-double RT0tri::LocalDiv(const Point &z, int i) const {
-    return 2.0;
-}
-
-VectorField RT0quad::LocalVector(const Point &z, int i) const {
-    switch (i) {
-        case 0:return VectorField(0.0, z[1] - 1.0);
-        case 1:return VectorField(z[0], 0.0);
-        case 2:return VectorField(0.0, z[1]);
-        case 3:return VectorField(z[0] - 1.0, 0.0);
-    }
-    Exit("Not implemented");
-    return 0.0;
-}
-
-double RT0quad::LocalDiv(const Point &z, int i) const {
-    return 1.0;
-}
-
-VectorField RT0tet::LocalVector(const Point &z, int i) const {
-    return zero;
-}
-
-double RT0tet::LocalDiv(const Point &z, int i) const {
-    return 0.0;
-}
-
-VectorField RT0hex::LocalVector(const Point &z, int i) const {
-    switch (i) {
-        case 0:return VectorField(0.0, 0.0, z[2] - 1.0);
-        case 1:return VectorField(0.0, z[1] - 1.0, 0.0);
-        case 2:return VectorField(z[0], 0.0, 0.0);
-        case 3:return VectorField(0.0, z[1], 0.0);
-        case 4:return VectorField(z[0] - 1.0, 0.0, 0.0);
-        case 5:return VectorField(0.0, 0.0, z[2]);
-    }
-    Exit("Not implemented");
-    return 0.0;
-}
-
 class P1rotquad: public Shape {
 public:
     const char *Name() const { return "P1rotquad"; }
@@ -1865,6 +1752,10 @@ Shape *Get_Shape(const string &name, const Quadrature &quadrature,
         return new P2curl_tet(quadrature, faceQuadrature);
     if (name == "CR_P1tri")
         return new CR_P1tri(quadrature, faceQuadrature);
+
+    // Raviart–Thomas
+    if (name == "RTOInt")
+        return new RT0Int(quadrature, faceQuadrature);
     if (name == "RT0tri")
         return new RT0tri(quadrature, faceQuadrature);
     if (name == "RT0quad")