From 2938548f8876e1a737875e539f16dec4a5e218f5 Mon Sep 17 00:00:00 2001
From: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date: Mon, 31 May 2010 12:38:28 +0000
Subject: [PATCH] dg : implicit 3d benchmark + order 6 to 10 tet polynomial
basis
---
Common/GmshDefines.h | 7 ++++++-
Geo/MTetrahedron.cpp | 5 +++++
Numeric/polynomialBasis.cpp | 30 ++++++++++++++++++++++++++++++
Solver/linearSystemPETSc.cpp | 2 +-
4 files changed, 42 insertions(+), 2 deletions(-)
diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index 6b0e1e201a..4859b4ef3c 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -126,8 +126,13 @@
#define MSH_TRI_B 68
#define MSH_POLYG_B 69
#define MSH_LIN_C 70
+#define MSH_TET_84 71
+#define MSH_TET_120 72
+#define MSH_TET_165 73
+#define MSH_TET_220 74
+#define MSH_TET_286 75
-#define MSH_NUM_TYPE 70
+#define MSH_NUM_TYPE 71
// Geometric entities
#define ENT_NONE 0
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index ea523860ec..ca1cc08c55 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -126,6 +126,11 @@ const polynomialBasis* MTetrahedron::getFunctionSpace(int o) const
case 3: return &polynomialBases::find(MSH_TET_20);
case 4: return &polynomialBases::find(MSH_TET_35);
case 5: return &polynomialBases::find(MSH_TET_56);
+ case 6: return &polynomialBases::find(MSH_TET_84);
+ case 7: return &polynomialBases::find(MSH_TET_120);
+ case 8: return &polynomialBases::find(MSH_TET_165);
+ case 9: return &polynomialBases::find(MSH_TET_220);
+ case 10: return &polynomialBases::find(MSH_TET_286);
default: Msg::Error("Order %d tetrahedron function space not implemented", order);
}
}
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 51c1cdee03..39b93f36d2 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -1114,6 +1114,36 @@ const polynomialBasis &polynomialBases::find(int tag)
F.points = gmshGeneratePointsTetrahedron(5, false);
generate3dFaceClosure(F.closures, 5);
break;
+ case MSH_TET_84 :
+ F.numFaces = 4;
+ F.monomials = generatePascalTetrahedron(6);
+ F.points = gmshGeneratePointsTetrahedron(6, false);
+ generate3dFaceClosure(F.closures, 6);
+ break;
+ case MSH_TET_120 :
+ F.numFaces = 4;
+ F.monomials = generatePascalTetrahedron(7);
+ F.points = gmshGeneratePointsTetrahedron(7, false);
+ generate3dFaceClosure(F.closures, 7);
+ break;
+ case MSH_TET_165 :
+ F.numFaces = 4;
+ F.monomials = generatePascalTetrahedron(8);
+ F.points = gmshGeneratePointsTetrahedron(8, false);
+ generate3dFaceClosure(F.closures, 8);
+ break;
+ case MSH_TET_220 :
+ F.numFaces = 4;
+ F.monomials = generatePascalTetrahedron(9);
+ F.points = gmshGeneratePointsTetrahedron(9, false);
+ generate3dFaceClosure(F.closures, 9);
+ break;
+ case MSH_TET_286 :
+ F.numFaces = 4;
+ F.monomials = generatePascalTetrahedron(10);
+ F.points = gmshGeneratePointsTetrahedron(10, false);
+ generate3dFaceClosure(F.closures, 10);
+ break;
case MSH_QUA_4 :
F.numFaces = 4;
F.monomials = generatePascalQuad(1);
diff --git a/Solver/linearSystemPETSc.cpp b/Solver/linearSystemPETSc.cpp
index 82805bf9a0..a2fdcdb777 100644
--- a/Solver/linearSystemPETSc.cpp
+++ b/Solver/linearSystemPETSc.cpp
@@ -83,7 +83,7 @@ void linearSystemPETSc<fullMatrix<PetscScalar> >::allocate(int nbRows)
// override the default options with the ones from the option
// database (if any)
_try(MatSetFromOptions(_a));
- _try(MatSeqBAIJSetPreallocation(_a, _blockSize, 4, NULL)); //todo preAllocate off-diagonal part
+ _try(MatSeqBAIJSetPreallocation(_a, _blockSize, 5, NULL)); //todo preAllocate off-diagonal part
//_try(MatMPIBAIJSetPreallocation(_a, _blockSize, 4, NULL, 0, NULL)); //todo preAllocate off-diagonal part
_try(VecCreate(PETSC_COMM_WORLD, &_x));
_try(VecSetSizes(_x, PETSC_DECIDE, nbRows * _blockSize));
--
GitLab