diff --git a/Geo/MElement.h b/Geo/MElement.h
index 3a6e92723d815b2c05cb91ba0edc24468b60dfea..5e5af4fc6410ac15caf566111d534b55c7894713 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -237,7 +237,7 @@ class MElement
// integration routines
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
{
- Msg::Error("No integration points defined for this type of element");
+ Msg::Error("No integration points defined for this type of element: %d",this->getType());
}
// IO routines
diff --git a/Geo/MPrism.cpp b/Geo/MPrism.cpp
index f4f30c5c75fdbd10237f7ae0d89e810f8c59bd33..dbe7e9b8ac2307fd7b1211a1ab59ad88c4f97c64 100644
--- a/Geo/MPrism.cpp
+++ b/Geo/MPrism.cpp
@@ -20,3 +20,91 @@ int MPrism::getVolumeSign()
mat[2][2] = _v[3]->z() - _v[0]->z();
return sign(det3x3(mat));
}
+
+void MPrism::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
+{
+ *npts = getNGQPriPts(pOrder);
+ *pts = getGQPriPts(pOrder);
+}
+
+const polynomialBasis* MPrism::getFunctionSpace(int o) const
+{
+ int order = (o == -1) ? getPolynomialOrder() : o;
+
+ int nv = getNumVolumeVertices();
+
+ if ((nv == 0) && (o == -1)) {
+ switch (order) {
+ case 1: return &polynomialBases::find(MSH_PRI_6);
+ case 2: return &polynomialBases::find(MSH_PRI_18);
+ default: Msg::Error("Order %d tetrahedron function space not implemented", order);
+ }
+ }
+ else {
+ switch (order) {
+ case 1: return &polynomialBases::find(MSH_PRI_6);
+ case 2: return &polynomialBases::find(MSH_PRI_18);
+ default: Msg::Error("Order %d tetrahedron function space not implemented", order);
+ }
+ }
+ return 0;
+}
+
+void MPrism::getFaceInfo(const MFace &face, int &ithFace, int &sign, int &rot) const
+{
+ for (ithFace = 0; ithFace < 5; ithFace++){
+ MVertex *v0 = _v[faces_prism(ithFace, 0)];
+ MVertex *v1 = _v[faces_prism(ithFace, 1)];
+ MVertex *v2 = _v[faces_prism(ithFace, 2)];
+
+ if (face.getNumVertices()==3) {
+ if (v0 == face.getVertex(0) && v1 == face.getVertex(1) && v2 == face.getVertex(2)){
+ sign = 1; rot = 0; return;
+ }
+ if (v0 == face.getVertex(1) && v1 == face.getVertex(2) && v2 == face.getVertex(0)){
+ sign = 1; rot = 1; return;
+ }
+ if (v0 == face.getVertex(2) && v1 == face.getVertex(0) && v2 == face.getVertex(1)){
+ sign = 1; rot = 2; return;
+ }
+ if (v0 == face.getVertex(0) && v1 == face.getVertex(2) && v2 == face.getVertex(1)){
+ sign = -1; rot = 0; return;
+ }
+ if (v0 == face.getVertex(1) && v1 == face.getVertex(0) && v2 == face.getVertex(2)){
+ sign = -1; rot = 1; return;
+ }
+ if (v0 == face.getVertex(2) && v1 == face.getVertex(1) && v2 == face.getVertex(0)){
+ sign = -1; rot = 2; return;
+ }
+ }
+ else {
+ MVertex *v3 = _v[faces_prism(ithFace, 3)];
+ if (v0 == face.getVertex(0) && v1 == face.getVertex(1) && v2 == face.getVertex(2) && v3 == face.getVertex(3)){
+ sign = 1; rot = 0; return;
+ }
+ if (v0 == face.getVertex(1) && v1 == face.getVertex(2) && v2 == face.getVertex(3) && v3 == face.getVertex(0)){
+ sign = 1; rot = 1; return;
+ }
+ if (v0 == face.getVertex(2) && v1 == face.getVertex(3) && v2 == face.getVertex(0) && v3 == face.getVertex(1)){
+ sign = 1; rot = 2; return;
+ }
+ if (v0 == face.getVertex(3) && v1 == face.getVertex(0) && v2 == face.getVertex(1) && v3 == face.getVertex(2)){
+ sign = 1; rot = 3; return;
+ }
+ if (v0 == face.getVertex(0) && v1 == face.getVertex(3) && v2 == face.getVertex(2) && v3 == face.getVertex(1)){
+ sign = -1; rot = 0; return;
+ }
+ if (v0 == face.getVertex(1) && v1 == face.getVertex(0) && v2 == face.getVertex(3) && v3 == face.getVertex(2)){
+ sign = -1; rot = 1; return;
+ }
+ if (v0 == face.getVertex(2) && v1 == face.getVertex(1) && v2 == face.getVertex(0) && v3 == face.getVertex(3)){
+ sign = -1; rot = 2; return;
+ }
+ if (v0 == face.getVertex(3) && v1 == face.getVertex(2) && v2 == face.getVertex(1) && v3 == face.getVertex(0)){
+ sign = -1; rot = 3; return;
+ }
+ }
+
+ }
+ Msg::Error("Could not get face information for prism %d", getNum());
+}
diff --git a/Geo/MPrism.h b/Geo/MPrism.h
index 01750011db3329d25c8e937cd8309d6fda9006e9..3d1220f472eeda675f435e2d3b853582edc44958 100644
--- a/Geo/MPrism.h
+++ b/Geo/MPrism.h
@@ -86,6 +86,7 @@ class MPrism : public MElement {
_getEdgeVertices(num, v);
}
virtual int getNumFaces(){ return 5; }
+ virtual void getFaceInfo(const MFace & face, int &ithFace, int &sign, int &rot) const;
virtual MFace getFace(int num)
{
if(num < 2)
@@ -128,6 +129,7 @@ class MPrism : public MElement {
tmp = _v[0]; _v[0] = _v[1]; _v[1] = tmp;
tmp = _v[3]; _v[3] = _v[4]; _v[4] = tmp;
}
+ virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
virtual int getVolumeSign();
virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
{
@@ -167,6 +169,7 @@ class MPrism : public MElement {
return false;
return true;
}
+ virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
private:
int edges_prism(const int edge, const int vert) const
{
diff --git a/Numeric/CMakeLists.txt b/Numeric/CMakeLists.txt
index f3ffff1f5e9822714f3e3a664de9838f25f1e648..77f7f0bc4dbec31d33ad272cab58cd18aa248df1 100644
--- a/Numeric/CMakeLists.txt
+++ b/Numeric/CMakeLists.txt
@@ -12,6 +12,7 @@ set(SRC
GaussQuadratureQuad.cpp
GaussQuadratureTet.cpp
GaussQuadratureHex.cpp
+ GaussQuadraturePri.cpp
GaussLegendreSimplex.cpp
robustPredicates.cpp
mathEvaluator.cpp
diff --git a/Numeric/Gauss.h b/Numeric/Gauss.h
index 6958223b65467d1fc3e65d74f7ed6cbfad64676c..d3f8791b59608c846a6054a7640887ed00008190 100644
--- a/Numeric/Gauss.h
+++ b/Numeric/Gauss.h
@@ -14,6 +14,7 @@ struct IntPt{
int GaussLegendreTri(int n1, int n2, IntPt *pts);
int GaussLegendreTet(int n1, int n2, int n3, IntPt *pts);
int GaussLegendreHex(int n1, int n2, int n3, IntPt *pts);
+int GaussLegendrePri(int n1, int n2, int n3, IntPt *pts);
int getNGQLPts (int order);
IntPt *getGQLPts (int order);
@@ -27,6 +28,9 @@ IntPt *getGQQPts(int order);
int getNGQTetPts(int order);
IntPt *getGQTetPts(int order);
+int getNGQPriPts(int order);
+IntPt *getGQPriPts(int order);
+
int getNGQHPts(int order);
IntPt *getGQHPts(int order);
diff --git a/Numeric/GaussQuadraturePri.cpp b/Numeric/GaussQuadraturePri.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3996d2954369d132fa98ac67479d50356e867f86
--- /dev/null
+++ b/Numeric/GaussQuadraturePri.cpp
@@ -0,0 +1,58 @@
+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#include "Gauss.h"
+#include "GaussLegendre1D.h"
+#include "stdio.h"
+
+IntPt *getGQPriPts(int order);
+int getNGQPriPts(int order);
+
+IntPt * GQP[20] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
+
+IntPt *getGQPriPts(int order)
+{
+ int nLin = (order+3)/2;
+ int nTri = getNGQTPts(order);
+ int n = nLin*nTri;
+ int index = n;
+ if(!GQP[index])
+ {
+ double *linPt,*linWt;
+ IntPt *triPts = getGQTPts(order);
+// printf("o = %d n = %d nT = %d nL = %d i = %d\n",order,n,nTri,nLin,index);
+ gmshGaussLegendre1D(nLin,&linPt,&linWt);
+ GQP[index] = new IntPt[n];
+ int l = 0;
+ for(int i=0; i < nTri; i++) {
+ for(int j=0; j < nLin; j++) {
+ GQP[index][l].pt[0] = triPts[i].pt[0];
+ GQP[index][l].pt[1] = triPts[i].pt[1];
+ GQP[index][l].pt[2] = linPt[j];
+ GQP[index][l++].weight = triPts[i].weight*linWt[j];
+// printf ("%d: %f %f %f %f\n",l-1,triPts[i].pt[0],triPts[i].pt[1],linPt[j],triPts[i].weight*linWt[j]);
+ }
+ }
+ }
+ return GQP[index];
+}
+
+int getNGQPriPts(int order)
+{
+ int nLin = (order+3)/2;
+ int nTri = getNGQTPts(order);
+ return nLin*nTri;
+
+// if(order == 3)return 8;
+// if(order == 2)return 8;
+// if(order < 2)
+// return GQPnPt[order];
+// return ((order+3)/2)*((order+3)/2)*((order+3)/2);
+}
+
+
+
+
+
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 2ddea17fc75fcea4524db2147ee4a17405131c36..537a4b2f43215dc911952660acf130dc9952336e 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -233,6 +233,32 @@ static fullMatrix<double> generatePascalTetrahedron(int order)
return monomials;
}
+// generate Pascal prism
+
+static fullMatrix<double> generatePascalPrism(int order)
+{
+ int nbMonomials = (order + 1) * (order + 1) * (order + 2) / 2;
+
+ fullMatrix<double> monomials(nbMonomials, 3);
+
+ int index = 0;
+ fullMatrix<double> lineMonoms = generate1DMonomials(order);
+ fullMatrix<double> triMonoms = generatePascalTriangle(order);
+// printf("nb: %d nT: %d %d nL: %d %d\n",nbMonomials,triMonoms.size1(),(order+1)*(order+2)/2,lineMonoms.size1(),order+1);
+ for (int j = 0; j < lineMonoms.size1(); j++) {
+ for (int i = 0; i < triMonoms.size1(); i++) {
+ monomials(index,0) = triMonoms(i,0);
+ monomials(index,1) = triMonoms(i,1);
+ monomials(index,2) = lineMonoms(j,0);
+ index ++;
+// printf("%d: %f %f %f\n",index,triMonoms(i,0),triMonoms(i,1),lineMonoms(j,0));
+ }
+ }
+// monomials.print();
+ return monomials;
+}
+
+
static int nbdoftriangle(int order) { return (order + 1) * (order + 2) / 2; }
//static int nbdoftriangleserendip(int order) { return 3 * order; }
@@ -597,6 +623,31 @@ static fullMatrix<double> gmshGeneratePointsTriangle(int order, bool serendip)
return point;
}
+static fullMatrix<double> gmshGeneratePointsPrism(int order, bool serendip)
+{
+ int nbPoints = (order + 1)*(order + 1)*(order + 2)/2;
+ fullMatrix<double> point(nbPoints, 3);
+
+ double overOrder = (order == 0 ? 1. : 1. / order);
+ int index = 0;
+ fullMatrix<double> triPoints = gmshGeneratePointsTriangle(order,false);
+ fullMatrix<double> linePoints = generate1DPoints(order);
+// printf("nb: %d nT: %d %d nL: %d %d\n",nbPoints,triPoints.size1(),(order+1)*(order+2)/2,linePoints.size1(),order+1);
+ for (int j = 0; j < linePoints.size1(); j++) {
+ for (int i = 0; i < triPoints.size1(); i++) {
+ point(index,0) = triPoints(i,0);
+ point(index,1) = triPoints(i,1);
+ point(index,2) = linePoints(j,0);
+ index ++;
+// printf("%d: %f %f %f\n",index,triPoints(i,0),triPoints(i,1),linePoints(j,0));
+ }
+ }
+// point.print();
+
+ point.scale(overOrder);
+ return point;
+}
+
static fullMatrix<double> gmshGeneratePointsQuad(int order, bool serendip)
{
int nbPoints = serendip ? order*4 : (order+1)*(order+1);
@@ -744,6 +795,50 @@ static void generate3dFaceClosure(polynomialBasis::clCont &closure, int order)
}
}
+static void getFaceClosurePrism(int iFace, int iSign, int iRotate, std::vector<int> &closure,
+ int order)
+{
+ if (order > 1)
+ Msg::Error("FaceClosure not implemented for prisms of order %d",order);
+ bool isTriangle = iFace<2;
+ int nNodes = isTriangle ? (order+1)*(order+2)/2 : (order+1)*(order+1);
+ closure.clear();
+ closure.resize(nNodes);
+ switch (order){
+ case 0:
+ closure[0] = 0;
+ break;
+ default:
+// int face[4][3] = {{-3, -2, -1}, {1, -6, 4}, {-4, 5, 3}, {6, 2, -5}};
+ int order1node[5][4] = {{0, 2, 1, -1}, {3, 4, 5, -1}, {0, 1, 4, 3}, {0, 3, 5, 2}, {1, 2, 5, 4}};
+ int nVertex = isTriangle ? 3 : 4;
+ for (int i = 0; i < nVertex; ++i){
+ int k;
+ k = (nVertex + (iSign * i) + iRotate) % nVertex; //- iSign * iRotate
+ closure[i] = order1node[iFace][k];
+ }
+ break;
+ }
+}
+
+static void generate3dFaceClosurePrism(polynomialBasis::clCont &closure, int order)
+{
+
+ closure.clear();
+ for (int iRotate = 0; iRotate < 4; iRotate++){
+ for (int iSign = 1; iSign >= -1; iSign -= 2){
+ for (int iFace = 0; iFace < 5; iFace++){
+ std::vector<int> closure_face;
+ getFaceClosurePrism(iFace, iSign, iRotate, closure_face, order);
+ closure.push_back(closure_face);
+// for(int i=0;i < closure_face.size(); i++) { printf("%d ",closure_face.at(i)); }
+// printf("\n");
+ }
+ }
+ }
+}
+
+
static void generate2dEdgeClosure(polynomialBasis::clCont &closure, int order, int nNod = 3)
{
closure.clear();
@@ -776,7 +871,8 @@ const polynomialBasis &polynomialBases::find(int tag)
std::map<int, polynomialBasis>::const_iterator it = fs.find(tag);
if (it != fs.end()) return it->second;
polynomialBasis F;
-
+ F.numFaces = -1;
+
switch (tag){
case MSH_PNT:
F.monomials = generate1DMonomials(0);
@@ -848,36 +944,43 @@ const polynomialBasis &polynomialBases::find(int tag)
generate2dEdgeClosure(F.edgeClosure, 5);
break;
case MSH_TET_4 :
+ F.numFaces = 4;
F.monomials = generatePascalTetrahedron(1);
F.points = gmshGeneratePointsTetrahedron(1, false);
generate3dFaceClosure(F.faceClosure, 1);
break;
case MSH_TET_10 :
+ F.numFaces = 4;
F.monomials = generatePascalTetrahedron(2);
F.points = gmshGeneratePointsTetrahedron(2, false);
generate3dFaceClosure(F.faceClosure, 2);
break;
case MSH_TET_20 :
+ F.numFaces = 4;
F.monomials = generatePascalTetrahedron(3);
F.points = gmshGeneratePointsTetrahedron(3, false);
generate3dFaceClosure(F.faceClosure, 3);
break;
case MSH_TET_35 :
+ F.numFaces = 4;
F.monomials = generatePascalTetrahedron(4);
F.points = gmshGeneratePointsTetrahedron(4, false);
generate3dFaceClosure(F.faceClosure, 4);
break;
case MSH_TET_34 :
+ F.numFaces = 4;
F.monomials = generatePascalSerendipityTetrahedron(4);
F.points = gmshGeneratePointsTetrahedron(4, true);
generate3dFaceClosure(F.faceClosure, 4);
break;
case MSH_TET_52 :
+ F.numFaces = 4;
F.monomials = generatePascalSerendipityTetrahedron(5);
F.points = gmshGeneratePointsTetrahedron(5, true);
generate3dFaceClosure(F.faceClosure, 5);
break;
case MSH_TET_56 :
+ F.numFaces = 4;
F.monomials = generatePascalTetrahedron(5);
F.points = gmshGeneratePointsTetrahedron(5, false);
generate3dFaceClosure(F.faceClosure, 5);
@@ -927,6 +1030,19 @@ const polynomialBasis &polynomialBases::find(int tag)
F.points = gmshGeneratePointsQuad(5,true);
generate2dEdgeClosure(F.edgeClosure, 5, 4);
break;
+ case MSH_PRI_6 : // first order
+ F.numFaces = 5;
+ F.monomials = generatePascalPrism(1);
+ F.points = gmshGeneratePointsPrism(1, false);
+ generate3dFaceClosurePrism(F.faceClosure, 1);
+ break;
+ case MSH_PRI_18 : // second order
+ F.numFaces = 5;
+ F.monomials = generatePascalPrism(2);
+ F.points = gmshGeneratePointsPrism(2, false);
+ generate3dFaceClosurePrism(F.faceClosure, 2);
+ break;
+
default :
Msg::Error("Unknown function space %d: reverting to TET_4", tag);
F.monomials = generatePascalTetrahedron(1);
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index f0705daf30f73b9b2cd10153ec00f8baa57dfeb9..ef0c624c992efff6d3463604d4397f787e77f772 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -23,10 +23,13 @@ class polynomialBasis
fullMatrix<double> points;
fullMatrix<double> monomials;
fullMatrix<double> coefficients;
+ int numFaces;
// for a given face/edge, with both a sign and a rotation,
// give an ordered list of nodes on this face/edge
inline int getFaceClosureId(int iFace, int iSign, int iRot) const {
- return iFace + 4 * (iSign == 1 ? 0 : 1) + 8 * iRot;
+// return iFace + 4 * (iSign == 1 ? 0 : 1) + 8 * iRot; // only for tetrahedra
+// return iFace + 5*(iSign == 1 ? 0 : 1) + 10*iRot; // only for prisms
+ return iFace + numFaces*(iSign == 1 ? 0 : 1) + 2*numFaces*iRot; // both tetrahedra and prisms
}
inline const std::vector<int> &getFaceClosure(int id) const
{
diff --git a/Solver/TESTCASES/Advection1D.lua b/Solver/TESTCASES/Advection1D.lua
index f8994e7e8bde69ccb239cad44b668cf6e282bc11..948af7171a8070882aa7f6ba9f924103690e35f7 100644
--- a/Solver/TESTCASES/Advection1D.lua
+++ b/Solver/TESTCASES/Advection1D.lua
@@ -16,7 +16,7 @@ v:set(2,0,0)
nu=fullMatrix(1,1);
nu:set(0,0,0)
-law = dgConservationLawAdvection(functionConstant(v):getName(), '')
+law = dgConservationLawAdvectionDiffusion(functionConstant(v):getName(), '')
--FunctionConstant(nu):getName())
dg:setConservationLaw(law)
diff --git a/Solver/TESTCASES/AdvectionDiffusion.lua b/Solver/TESTCASES/AdvectionDiffusion.lua
index 56a34529dc1daa5bc7af4b3cd777006d16dd3c13..13b691ddc757c4359d269d27cea0d262e210c657 100644
--- a/Solver/TESTCASES/AdvectionDiffusion.lua
+++ b/Solver/TESTCASES/AdvectionDiffusion.lua
@@ -24,7 +24,7 @@ v:set(2,0,0)
nu=fullMatrix(1,1);
nu:set(0,0,0.001)
-law = dgConservationLawAdvection(functionConstant(v):getName(),functionConstant(nu):getName())
+law = dgConservationLawAdvectionDiffusion(functionConstant(v):getName(),functionConstant(nu):getName())
dg:setConservationLaw(law)
-- boundary condition
diff --git a/Solver/TESTCASES/edge.msh b/Solver/TESTCASES/edge.msh
index 5ee49d3930ffa1973e9913f8c37f65abe2d39865..b77e58a920be6af1bbca0e3bc747580e9c5107a1 100644
--- a/Solver/TESTCASES/edge.msh
+++ b/Solver/TESTCASES/edge.msh
@@ -1,5 +1,5 @@
$MeshFormat
-2.1 0 8
+2.2 0 8
$EndMeshFormat
$PhysicalNames
3
@@ -212,205 +212,205 @@ $Nodes
$EndNodes
$Elements
201
-1 15 3 1 1 0 1
-2 15 3 2 2 0 2
-3 1 3 3 1 0 1 3
-4 1 3 3 1 0 3 4
-5 1 3 3 1 0 4 5
-6 1 3 3 1 0 5 6
-7 1 3 3 1 0 6 7
-8 1 3 3 1 0 7 8
-9 1 3 3 1 0 8 9
-10 1 3 3 1 0 9 10
-11 1 3 3 1 0 10 11
-12 1 3 3 1 0 11 12
-13 1 3 3 1 0 12 13
-14 1 3 3 1 0 13 14
-15 1 3 3 1 0 14 15
-16 1 3 3 1 0 15 16
-17 1 3 3 1 0 16 17
-18 1 3 3 1 0 17 18
-19 1 3 3 1 0 18 19
-20 1 3 3 1 0 19 20
-21 1 3 3 1 0 20 21
-22 1 3 3 1 0 21 22
-23 1 3 3 1 0 22 23
-24 1 3 3 1 0 23 24
-25 1 3 3 1 0 24 25
-26 1 3 3 1 0 25 26
-27 1 3 3 1 0 26 27
-28 1 3 3 1 0 27 28
-29 1 3 3 1 0 28 29
-30 1 3 3 1 0 29 30
-31 1 3 3 1 0 30 31
-32 1 3 3 1 0 31 32
-33 1 3 3 1 0 32 33
-34 1 3 3 1 0 33 34
-35 1 3 3 1 0 34 35
-36 1 3 3 1 0 35 36
-37 1 3 3 1 0 36 37
-38 1 3 3 1 0 37 38
-39 1 3 3 1 0 38 39
-40 1 3 3 1 0 39 40
-41 1 3 3 1 0 40 41
-42 1 3 3 1 0 41 42
-43 1 3 3 1 0 42 43
-44 1 3 3 1 0 43 44
-45 1 3 3 1 0 44 45
-46 1 3 3 1 0 45 46
-47 1 3 3 1 0 46 47
-48 1 3 3 1 0 47 48
-49 1 3 3 1 0 48 49
-50 1 3 3 1 0 49 50
-51 1 3 3 1 0 50 51
-52 1 3 3 1 0 51 52
-53 1 3 3 1 0 52 53
-54 1 3 3 1 0 53 54
-55 1 3 3 1 0 54 55
-56 1 3 3 1 0 55 56
-57 1 3 3 1 0 56 57
-58 1 3 3 1 0 57 58
-59 1 3 3 1 0 58 59
-60 1 3 3 1 0 59 60
-61 1 3 3 1 0 60 61
-62 1 3 3 1 0 61 62
-63 1 3 3 1 0 62 63
-64 1 3 3 1 0 63 64
-65 1 3 3 1 0 64 65
-66 1 3 3 1 0 65 66
-67 1 3 3 1 0 66 67
-68 1 3 3 1 0 67 68
-69 1 3 3 1 0 68 69
-70 1 3 3 1 0 69 70
-71 1 3 3 1 0 70 71
-72 1 3 3 1 0 71 72
-73 1 3 3 1 0 72 73
-74 1 3 3 1 0 73 74
-75 1 3 3 1 0 74 75
-76 1 3 3 1 0 75 76
-77 1 3 3 1 0 76 77
-78 1 3 3 1 0 77 78
-79 1 3 3 1 0 78 79
-80 1 3 3 1 0 79 80
-81 1 3 3 1 0 80 81
-82 1 3 3 1 0 81 82
-83 1 3 3 1 0 82 83
-84 1 3 3 1 0 83 84
-85 1 3 3 1 0 84 85
-86 1 3 3 1 0 85 86
-87 1 3 3 1 0 86 87
-88 1 3 3 1 0 87 88
-89 1 3 3 1 0 88 89
-90 1 3 3 1 0 89 90
-91 1 3 3 1 0 90 91
-92 1 3 3 1 0 91 92
-93 1 3 3 1 0 92 93
-94 1 3 3 1 0 93 94
-95 1 3 3 1 0 94 95
-96 1 3 3 1 0 95 96
-97 1 3 3 1 0 96 97
-98 1 3 3 1 0 97 98
-99 1 3 3 1 0 98 99
-100 1 3 3 1 0 99 100
-101 1 3 3 1 0 100 101
-102 1 3 3 1 0 101 102
-103 1 3 3 1 0 102 103
-104 1 3 3 1 0 103 104
-105 1 3 3 1 0 104 105
-106 1 3 3 1 0 105 106
-107 1 3 3 1 0 106 107
-108 1 3 3 1 0 107 108
-109 1 3 3 1 0 108 109
-110 1 3 3 1 0 109 110
-111 1 3 3 1 0 110 111
-112 1 3 3 1 0 111 112
-113 1 3 3 1 0 112 113
-114 1 3 3 1 0 113 114
-115 1 3 3 1 0 114 115
-116 1 3 3 1 0 115 116
-117 1 3 3 1 0 116 117
-118 1 3 3 1 0 117 118
-119 1 3 3 1 0 118 119
-120 1 3 3 1 0 119 120
-121 1 3 3 1 0 120 121
-122 1 3 3 1 0 121 122
-123 1 3 3 1 0 122 123
-124 1 3 3 1 0 123 124
-125 1 3 3 1 0 124 125
-126 1 3 3 1 0 125 126
-127 1 3 3 1 0 126 127
-128 1 3 3 1 0 127 128
-129 1 3 3 1 0 128 129
-130 1 3 3 1 0 129 130
-131 1 3 3 1 0 130 131
-132 1 3 3 1 0 131 132
-133 1 3 3 1 0 132 133
-134 1 3 3 1 0 133 134
-135 1 3 3 1 0 134 135
-136 1 3 3 1 0 135 136
-137 1 3 3 1 0 136 137
-138 1 3 3 1 0 137 138
-139 1 3 3 1 0 138 139
-140 1 3 3 1 0 139 140
-141 1 3 3 1 0 140 141
-142 1 3 3 1 0 141 142
-143 1 3 3 1 0 142 143
-144 1 3 3 1 0 143 144
-145 1 3 3 1 0 144 145
-146 1 3 3 1 0 145 146
-147 1 3 3 1 0 146 147
-148 1 3 3 1 0 147 148
-149 1 3 3 1 0 148 149
-150 1 3 3 1 0 149 150
-151 1 3 3 1 0 150 151
-152 1 3 3 1 0 151 152
-153 1 3 3 1 0 152 153
-154 1 3 3 1 0 153 154
-155 1 3 3 1 0 154 155
-156 1 3 3 1 0 155 156
-157 1 3 3 1 0 156 157
-158 1 3 3 1 0 157 158
-159 1 3 3 1 0 158 159
-160 1 3 3 1 0 159 160
-161 1 3 3 1 0 160 161
-162 1 3 3 1 0 161 162
-163 1 3 3 1 0 162 163
-164 1 3 3 1 0 163 164
-165 1 3 3 1 0 164 165
-166 1 3 3 1 0 165 166
-167 1 3 3 1 0 166 167
-168 1 3 3 1 0 167 168
-169 1 3 3 1 0 168 169
-170 1 3 3 1 0 169 170
-171 1 3 3 1 0 170 171
-172 1 3 3 1 0 171 172
-173 1 3 3 1 0 172 173
-174 1 3 3 1 0 173 174
-175 1 3 3 1 0 174 175
-176 1 3 3 1 0 175 176
-177 1 3 3 1 0 176 177
-178 1 3 3 1 0 177 178
-179 1 3 3 1 0 178 179
-180 1 3 3 1 0 179 180
-181 1 3 3 1 0 180 181
-182 1 3 3 1 0 181 182
-183 1 3 3 1 0 182 183
-184 1 3 3 1 0 183 184
-185 1 3 3 1 0 184 185
-186 1 3 3 1 0 185 186
-187 1 3 3 1 0 186 187
-188 1 3 3 1 0 187 188
-189 1 3 3 1 0 188 189
-190 1 3 3 1 0 189 190
-191 1 3 3 1 0 190 191
-192 1 3 3 1 0 191 192
-193 1 3 3 1 0 192 193
-194 1 3 3 1 0 193 194
-195 1 3 3 1 0 194 195
-196 1 3 3 1 0 195 196
-197 1 3 3 1 0 196 197
-198 1 3 3 1 0 197 198
-199 1 3 3 1 0 198 199
-200 1 3 3 1 0 199 200
-201 1 3 3 1 0 200 2
+1 15 2 1 1 1
+2 15 2 2 2 2
+3 1 2 3 1 1 3
+4 1 2 3 1 3 4
+5 1 2 3 1 4 5
+6 1 2 3 1 5 6
+7 1 2 3 1 6 7
+8 1 2 3 1 7 8
+9 1 2 3 1 8 9
+10 1 2 3 1 9 10
+11 1 2 3 1 10 11
+12 1 2 3 1 11 12
+13 1 2 3 1 12 13
+14 1 2 3 1 13 14
+15 1 2 3 1 14 15
+16 1 2 3 1 15 16
+17 1 2 3 1 16 17
+18 1 2 3 1 17 18
+19 1 2 3 1 18 19
+20 1 2 3 1 19 20
+21 1 2 3 1 20 21
+22 1 2 3 1 21 22
+23 1 2 3 1 22 23
+24 1 2 3 1 23 24
+25 1 2 3 1 24 25
+26 1 2 3 1 25 26
+27 1 2 3 1 26 27
+28 1 2 3 1 27 28
+29 1 2 3 1 28 29
+30 1 2 3 1 29 30
+31 1 2 3 1 30 31
+32 1 2 3 1 31 32
+33 1 2 3 1 32 33
+34 1 2 3 1 33 34
+35 1 2 3 1 34 35
+36 1 2 3 1 35 36
+37 1 2 3 1 36 37
+38 1 2 3 1 37 38
+39 1 2 3 1 38 39
+40 1 2 3 1 39 40
+41 1 2 3 1 40 41
+42 1 2 3 1 41 42
+43 1 2 3 1 42 43
+44 1 2 3 1 43 44
+45 1 2 3 1 44 45
+46 1 2 3 1 45 46
+47 1 2 3 1 46 47
+48 1 2 3 1 47 48
+49 1 2 3 1 48 49
+50 1 2 3 1 49 50
+51 1 2 3 1 50 51
+52 1 2 3 1 51 52
+53 1 2 3 1 52 53
+54 1 2 3 1 53 54
+55 1 2 3 1 54 55
+56 1 2 3 1 55 56
+57 1 2 3 1 56 57
+58 1 2 3 1 57 58
+59 1 2 3 1 58 59
+60 1 2 3 1 59 60
+61 1 2 3 1 60 61
+62 1 2 3 1 61 62
+63 1 2 3 1 62 63
+64 1 2 3 1 63 64
+65 1 2 3 1 64 65
+66 1 2 3 1 65 66
+67 1 2 3 1 66 67
+68 1 2 3 1 67 68
+69 1 2 3 1 68 69
+70 1 2 3 1 69 70
+71 1 2 3 1 70 71
+72 1 2 3 1 71 72
+73 1 2 3 1 72 73
+74 1 2 3 1 73 74
+75 1 2 3 1 74 75
+76 1 2 3 1 75 76
+77 1 2 3 1 76 77
+78 1 2 3 1 77 78
+79 1 2 3 1 78 79
+80 1 2 3 1 79 80
+81 1 2 3 1 80 81
+82 1 2 3 1 81 82
+83 1 2 3 1 82 83
+84 1 2 3 1 83 84
+85 1 2 3 1 84 85
+86 1 2 3 1 85 86
+87 1 2 3 1 86 87
+88 1 2 3 1 87 88
+89 1 2 3 1 88 89
+90 1 2 3 1 89 90
+91 1 2 3 1 90 91
+92 1 2 3 1 91 92
+93 1 2 3 1 92 93
+94 1 2 3 1 93 94
+95 1 2 3 1 94 95
+96 1 2 3 1 95 96
+97 1 2 3 1 96 97
+98 1 2 3 1 97 98
+99 1 2 3 1 98 99
+100 1 2 3 1 99 100
+101 1 2 3 1 100 101
+102 1 2 3 1 101 102
+103 1 2 3 1 102 103
+104 1 2 3 1 103 104
+105 1 2 3 1 104 105
+106 1 2 3 1 105 106
+107 1 2 3 1 106 107
+108 1 2 3 1 107 108
+109 1 2 3 1 108 109
+110 1 2 3 1 109 110
+111 1 2 3 1 110 111
+112 1 2 3 1 111 112
+113 1 2 3 1 112 113
+114 1 2 3 1 113 114
+115 1 2 3 1 114 115
+116 1 2 3 1 115 116
+117 1 2 3 1 116 117
+118 1 2 3 1 117 118
+119 1 2 3 1 118 119
+120 1 2 3 1 119 120
+121 1 2 3 1 120 121
+122 1 2 3 1 121 122
+123 1 2 3 1 122 123
+124 1 2 3 1 123 124
+125 1 2 3 1 124 125
+126 1 2 3 1 125 126
+127 1 2 3 1 126 127
+128 1 2 3 1 127 128
+129 1 2 3 1 128 129
+130 1 2 3 1 129 130
+131 1 2 3 1 130 131
+132 1 2 3 1 131 132
+133 1 2 3 1 132 133
+134 1 2 3 1 133 134
+135 1 2 3 1 134 135
+136 1 2 3 1 135 136
+137 1 2 3 1 136 137
+138 1 2 3 1 137 138
+139 1 2 3 1 138 139
+140 1 2 3 1 139 140
+141 1 2 3 1 140 141
+142 1 2 3 1 141 142
+143 1 2 3 1 142 143
+144 1 2 3 1 143 144
+145 1 2 3 1 144 145
+146 1 2 3 1 145 146
+147 1 2 3 1 146 147
+148 1 2 3 1 147 148
+149 1 2 3 1 148 149
+150 1 2 3 1 149 150
+151 1 2 3 1 150 151
+152 1 2 3 1 151 152
+153 1 2 3 1 152 153
+154 1 2 3 1 153 154
+155 1 2 3 1 154 155
+156 1 2 3 1 155 156
+157 1 2 3 1 156 157
+158 1 2 3 1 157 158
+159 1 2 3 1 158 159
+160 1 2 3 1 159 160
+161 1 2 3 1 160 161
+162 1 2 3 1 161 162
+163 1 2 3 1 162 163
+164 1 2 3 1 163 164
+165 1 2 3 1 164 165
+166 1 2 3 1 165 166
+167 1 2 3 1 166 167
+168 1 2 3 1 167 168
+169 1 2 3 1 168 169
+170 1 2 3 1 169 170
+171 1 2 3 1 170 171
+172 1 2 3 1 171 172
+173 1 2 3 1 172 173
+174 1 2 3 1 173 174
+175 1 2 3 1 174 175
+176 1 2 3 1 175 176
+177 1 2 3 1 176 177
+178 1 2 3 1 177 178
+179 1 2 3 1 178 179
+180 1 2 3 1 179 180
+181 1 2 3 1 180 181
+182 1 2 3 1 181 182
+183 1 2 3 1 182 183
+184 1 2 3 1 183 184
+185 1 2 3 1 184 185
+186 1 2 3 1 185 186
+187 1 2 3 1 186 187
+188 1 2 3 1 187 188
+189 1 2 3 1 188 189
+190 1 2 3 1 189 190
+191 1 2 3 1 190 191
+192 1 2 3 1 191 192
+193 1 2 3 1 192 193
+194 1 2 3 1 193 194
+195 1 2 3 1 194 195
+196 1 2 3 1 195 196
+197 1 2 3 1 196 197
+198 1 2 3 1 197 198
+199 1 2 3 1 198 199
+200 1 2 3 1 199 200
+201 1 2 3 1 200 2
$EndElements
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index 837cd39486eb47031ac6d7ba8d27747520f70222..4f68cb58a30a4bf529421d97557e5fdf70c530d1 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -75,7 +75,7 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e,
_dPsiDx = new fullMatrix<double> ( _integration->size1()*3,nbPsi*e.size());
_elementVolume = new fullMatrix<double> (e.size(),1);
double g[256][3],f[256];
-
+
for (int i=0;i<_elements.size();i++){
MElement *e = _elements[i];
element_to_index[e] = i;
@@ -408,7 +408,7 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string bo
init(pOrder);
}
-dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, int pOrder):
+dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, int pOrder, int numVertices):
_groupLeft(elGroup),_groupRight(elGroup)
{
_fsLeft=_groupLeft.getElement(0)->getFunctionSpace(pOrder);
@@ -456,10 +456,12 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, int pOrder):
MElement &el = *elGroup.getElement(i);
for (int j=0; j<el.getNumFaces(); j++){
MFace face = el.getFace(j);
- if(faceMap.find(face) == faceMap.end()){
- faceMap[face] = i;
- }else{
- addFace(face,faceMap[face],i);
+ if (numVertices < 0 || face.getNumVertices() == numVertices) {
+ if(faceMap.find(face) == faceMap.end()){
+ faceMap[face] = i;
+ }else{
+ addFace(face,faceMap[face],i);
+ }
}
}
}
@@ -470,7 +472,7 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, int pOrder):
init(pOrder);
}
-dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup1, const dgGroupOfElements &elGroup2, int pOrder):
+dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup1, const dgGroupOfElements &elGroup2, int pOrder, int numVertices):
_groupLeft(elGroup1),_groupRight(elGroup2)
{
_fsLeft=_groupLeft.getElement(0)->getFunctionSpace(pOrder);
@@ -546,10 +548,12 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup1, const dgGroup
MElement &el = *elGroup1.getElement(i);
for (int j=0; j<el.getNumFaces(); j++){
MFace face = el.getFace(j);
- if(faceMap.find(face) == faceMap.end()){
- faceMap[face] = i;
- }else{
- faceMap.erase(face);
+ if (numVertices < 0 || face.getNumVertices() == numVertices) {
+ if(faceMap.find(face) == faceMap.end()){
+ faceMap[face] = i;
+ }else{
+ faceMap.erase(face);
+ }
}
}
}
@@ -703,6 +707,36 @@ void dgGroupCollection::buildGroupsOfElements(GModel *model, int dim, int order)
if(el->getPartition()==Msg::GetCommRank()+1 || el->getPartition()==0){
localElements[el->getType()].push_back(el);
nlocal++;
+ switch(dim) {
+ case 1: {
+ int interfaceType = 1; // number of vertices
+ if(_interfaceTypes.find(interfaceType) == _interfaceTypes.end()) {
+ _interfaceTypes.insert(interfaceType);
+// printf("Inserted interfaceType: %d el: %d dim: %d\n",interfaceType,el->getType(),dim);
+ }
+ break;
+ }
+ case 2: {
+ int interfaceType = 2; // number of vertices
+ if(_interfaceTypes.find(interfaceType) == _interfaceTypes.end()) {
+ _interfaceTypes.insert(interfaceType);
+// printf("Inserted interfaceType: %d el: %d dim: %d\n",interfaceType,el->getType(),dim);
+ }
+ break;
+ }
+ case 3: {
+ for(int iFace=0; iFace<el->getNumFaces(); iFace++) {
+ MFace face = el->getFace(iFace);
+ if(_interfaceTypes.find(face.getNumVertices()) == _interfaceTypes.end()) {
+ _interfaceTypes.insert(face.getNumVertices());
+// printf("Inserted interfaceType: %d el: %d dim: %d\n",face.getNumVertices(),el->getType(),dim);
+ }
+ }
+ break;
+ }
+ default :
+ throw;
+ }
}
}
}
@@ -730,7 +764,7 @@ void dgGroupCollection::buildGroupsOfInterfaces() {
std::map<const std::string,std::set<MVertex*> > boundaryVertices;
std::map<const std::string,std::set<MEdge, Less_Edge> > boundaryEdges;
- std::map<const std::string,std::set<MFace, Less_Face> > boundaryFaces;
+ std::map<const int,std::map<const std::string,std::set<MFace, Less_Face> > > boundaryFaces; // [elementType][bndString]
std::vector<std::map<int, std::vector<MElement *> > >ghostElements(Msg::GetCommSize()); // [partition][elementType]
int nghosts=0;
@@ -754,9 +788,11 @@ void dgGroupCollection::buildGroupsOfInterfaces() {
case 2:
boundaryEdges[physicalName].insert( element->getEdge(0) );
break;
- case 3:
- boundaryFaces[physicalName].insert( element->getFace(0));
- break;
+ case 3: {
+ MFace face = element->getFace(0);
+ boundaryFaces[element->getType()][physicalName].insert( face );
+ break;
+ }
default :
throw;
}
@@ -807,30 +843,36 @@ void dgGroupCollection::buildGroupsOfInterfaces() {
break;
}
case 3 : {
- std::map<const std::string, std::set<MFace, Less_Face> >::iterator mapIt;
- for(mapIt=boundaryFaces.begin(); mapIt!=boundaryFaces.end(); mapIt++) {
- dgGroupOfFaces *gof=new dgGroupOfFaces(*_elementGroups[i],mapIt->first,order,mapIt->second);
- if(gof->getNbElements())
- _boundaryGroups.push_back(gof);
- else
- delete gof;
+ std::map<const int,std::map<const std::string,std::set<MFace, Less_Face> > >::iterator elemTypeIt;
+ std::map<const std::string, std::set<MFace, Less_Face> >::iterator mapIt;
+ for(elemTypeIt=boundaryFaces.begin(); elemTypeIt!=boundaryFaces.end(); elemTypeIt++) {
+ for(mapIt=elemTypeIt->second.begin(); mapIt!=elemTypeIt->second.end(); mapIt++) {
+ dgGroupOfFaces *gof=new dgGroupOfFaces(*_elementGroups[i],mapIt->first,order,mapIt->second);
+ if(gof->getNbElements())
+ _boundaryGroups.push_back(gof);
+ else
+ delete gof;
+ }
}
break;
}
}
- // create a group of faces for every face that is common to elements of the same group
- dgGroupOfFaces *gof=new dgGroupOfFaces(*_elementGroups[i],order);
- if(gof->getNbElements())
- _faceGroups.push_back(gof);
- else
- delete gof;
- // create a groupe of d
- for (int j=i+1;j<_elementGroups.size();j++){
- dgGroupOfFaces *gof = new dgGroupOfFaces(*_elementGroups[i],*_elementGroups[j],order);
+ // create a group of faces for every face that is common to elements of the same group
+ // create separate groups for each face type
+ for(std::set<int>::iterator faceTypeIt=_interfaceTypes.begin(); faceTypeIt!=_interfaceTypes.end(); faceTypeIt++) {
+ dgGroupOfFaces *gof = new dgGroupOfFaces(*_elementGroups[i],order,*faceTypeIt);
if (gof->getNbElements())
_faceGroups.push_back(gof);
else
delete gof;
+// create a groupe of d
+ for (int j=i+1;j<_elementGroups.size();j++){
+ dgGroupOfFaces *gof = new dgGroupOfFaces(*_elementGroups[i],*_elementGroups[j],order,*faceTypeIt);
+ if (gof->getNbElements())
+ _faceGroups.push_back(gof);
+ else
+ delete gof;
+ }
}
}
//create ghost groups
@@ -846,14 +888,16 @@ void dgGroupCollection::buildGroupsOfInterfaces() {
//create face group for ghostGroups
for (int i=0; i<_ghostGroups.size(); i++){
for (int j=0;j<_elementGroups.size();j++){
- dgGroupOfFaces *gof = new dgGroupOfFaces(*_ghostGroups[i],*_elementGroups[j],order);
- if (gof->getNbElements())
- _faceGroups.push_back(gof);
- else
- delete gof;
+ for(std::set<int>::iterator faceTypeIt=_interfaceTypes.begin(); faceTypeIt!=_interfaceTypes.end(); faceTypeIt++) {
+ dgGroupOfFaces *gof = new dgGroupOfFaces(*_ghostGroups[i],*_elementGroups[j],order,*faceTypeIt);
+ if (gof->getNbElements())
+ _faceGroups.push_back(gof);
+ else
+ delete gof;
+ }
}
}
-
+ Msg::Info("%d groups of interfaces",_faceGroups.size());
// build the ghosts structure
int *nGhostElements = new int[Msg::GetCommSize()];
int *nParentElements = new int[Msg::GetCommSize()];
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 454403bbfca84816de6c25c0c0cd8d4982ad0c98..042a48c909e776d321bd5ea974ca498b819984b4 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -184,8 +184,8 @@ public:
inline fullMatrix<double> &getIntegrationOnElementRight(int iFace) { return _integrationPointsRight[_closuresIdRight[iFace]];}
inline fullMatrix<double> &getNormals () const {return *_normals;}
- dgGroupOfFaces (const dgGroupOfElements &elements,int pOrder);
- dgGroupOfFaces (const dgGroupOfElements &a, const dgGroupOfElements &b,int pOrder);
+ dgGroupOfFaces (const dgGroupOfElements &elements,int pOrder, int numVertices = -1);
+ dgGroupOfFaces (const dgGroupOfElements &a, const dgGroupOfElements &b,int pOrder, int numVertices = -1);
dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MVertex*> &boundaryVertices);
dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MEdge,Less_Edge> &boundaryEdges);
dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MFace,Less_Face> &boundaryFaces);
@@ -217,6 +217,8 @@ class dgGroupCollection {
std::vector<dgGroupOfFaces*> _faceGroups; //interface
std::vector<dgGroupOfFaces*> _boundaryGroups; //boundary
std::vector<dgGroupOfElements*> _ghostGroups; //ghost volume
+ // container of different face types (identified by the number of vertices)
+ std::set<int> _interfaceTypes;
//{group,id} of the elements to send to each partition for a scatter operation
std::vector< std::vector<std::pair<int,int> > >_elementsToSend;