From f413339c1af30d624f8d8b2230e210535fd4048d Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Wed, 18 May 2011 06:26:50 +0000
Subject: [PATCH] pp
---
Numeric/polynomialBasis.cpp | 78 +++++++++++++++----------------------
1 file changed, 31 insertions(+), 47 deletions(-)
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index d77b0644ee..3ff790f952 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -8,12 +8,12 @@
// Jonathan Lambrechts
//
+#include <stdlib.h>
#include "GmshDefines.h"
#include "GmshMessage.h"
#include "polynomialBasis.h"
#include "Gauss.h"
-#include "stdlib.h"
static void printClosure(polynomialBasis::clCont &fullClosure, std::vector<int> &closureRef,
polynomialBasis::clCont &closures)
{
@@ -191,8 +191,8 @@ static fullMatrix<double> generatePascalSerendipityTriangle(int order)
return monomials;
}
-const fullMatrix<double> &polynomialBasis::getGradientAtFaceIntegrationPoints(int integrationOrder,
- int closureId) const
+const fullMatrix<double> &polynomialBasis::
+getGradientAtFaceIntegrationPoints(int integrationOrder, int closureId) const
{
std::vector<fullMatrix<double> > &dfAtFace = _dfAtFace[integrationOrder];
if (dfAtFace.empty()) {
@@ -229,6 +229,7 @@ const fullMatrix<double> &polynomialBasis::getGradientAtFaceIntegrationPoints(in
}
// generate all monomials xi^m * eta^n with n and m <= order
+
static fullMatrix<double> generatePascalQuad(int order)
{
@@ -247,7 +248,6 @@ static fullMatrix<double> generatePascalQuad(int order)
return monomials;
}
-
/*
00 10 20 30 40 ⋯
01 11 21 31 41 ⋯
@@ -258,6 +258,7 @@ static fullMatrix<double> generatePascalQuad(int order)
*/
// generate all monomials xi^m * eta^n with n and m <= order
+
static fullMatrix<double> generatePascalHex(int order, bool serendip)
{
int siz = (order+1)*(order+1)*(order+1);
@@ -276,7 +277,6 @@ static fullMatrix<double> generatePascalHex(int order, bool serendip)
}
}
}
- // printf("siz %d %d\n",siz,index );
return monomials;
}
@@ -309,29 +309,6 @@ static fullMatrix<double> generatePascalQuadSerendip(int order)
return monomials;
}
-
-/*static fullMatrix<double> generatePascalQuadSerendip(int order)
-{
-
- fullMatrix<double> monomials( order*4, 2);
- int index = 0;
- for (int p = 0; p < order; p++) {
- monomials(p, 0) = p;
- monomials(p, 1) = 0;
-
- monomials(p+order, 0) = order;
- monomials(p+order, 1) = p;
-
- monomials(p+3*order, 0) = order-p;
- monomials(p+3*order, 1) = order;
-
- monomials(p+2*order, 0) = 0;
- monomials(p+2*order, 1) = order-p;
- }
- monomials.print();
- return monomials;
-}*/
-
// generate the monomials subspace of all monomials of order exactly == p
static fullMatrix<double> generateMonomialSubspace(int dim, int p)
@@ -447,12 +424,11 @@ static fullMatrix<double> generatePascalPrism(int order)
}
}
}
-// monomials.print("Pri monoms");
+
return monomials;
}
static int nbdoftriangle(int order) { return (order + 1) * (order + 2) / 2; }
-//static int nbdoftriangleserendip(int order) { return 3 * order; }
//KH : caveat : node coordinates are not yet coherent with node numbering associated
// to numbering of principal vertices of face !!!!
@@ -1055,8 +1031,8 @@ static void generate1dVertexClosure(polynomialBasis::clCont &closure, int order)
closure[1].type = MSH_PNT;
}
-static void generate1dVertexClosureFull(polynomialBasis::clCont &closure, std::vector<int> &closureRef,
- int order)
+static void generate1dVertexClosureFull(polynomialBasis::clCont &closure,
+ std::vector<int> &closureRef, int order)
{
closure.clear();
closure.resize(2);
@@ -1075,10 +1051,9 @@ static void generate1dVertexClosureFull(polynomialBasis::clCont &closure, std::v
closureRef[1] = 0;
}
-static void getFaceClosureTet(int iFace, int iSign, int iRotate, polynomialBasis::closure &closure,
- int order)
+static void getFaceClosureTet(int iFace, int iSign, int iRotate,
+ polynomialBasis::closure &closure, int order)
{
-
closure.clear();
closure.resize((order + 1) * (order + 2) / 2);
closure.type = polynomialBasis::getTag(TYPE_TRI, order, false);
@@ -1141,7 +1116,9 @@ static void getFaceClosureTet(int iFace, int iSign, int iRotate, polynomialBasis
break;
}
}
-static void generate2dEdgeClosureFull(polynomialBasis::clCont &closure, std::vector<int> &closureRef,
+
+static void generate2dEdgeClosureFull(polynomialBasis::clCont &closure,
+ std::vector<int> &closureRef,
int order, int nNod, bool serendip)
{
closure.clear();
@@ -1220,7 +1197,8 @@ static void generateFaceClosureTet(polynomialBasis::clCont &closure, int order)
}
}
-static void generateFaceClosureTetFull(polynomialBasis::clCont &closureFull, std::vector<int> &closureRef,
+static void generateFaceClosureTetFull(polynomialBasis::clCont &closureFull,
+ std::vector<int> &closureRef,
int order, bool serendip)
{
closureFull.clear();
@@ -1563,13 +1541,15 @@ const polynomialBasis *polynomialBases::find(int tag)
case TYPE_TRI :
F.numFaces = 3;
F.dimension = 2;
- F.monomials = F.serendip ? generatePascalSerendipityTriangle(F.order) : generatePascalTriangle(F.order);
+ F.monomials = F.serendip ? generatePascalSerendipityTriangle(F.order) :
+ generatePascalTriangle(F.order);
F.points = gmshGeneratePointsTriangle(F.order, F.serendip);
if (F.order == 0) {
generateClosureOrder0(F.closures, 6);
generateClosureOrder0(F.fullClosures, 6);
F.closureRef.resize(6, 0);
- } else {
+ }
+ else {
generate2dEdgeClosure(F.closures, F.order);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, F.order, 3, F.serendip);
}
@@ -1577,13 +1557,15 @@ const polynomialBasis *polynomialBases::find(int tag)
case TYPE_QUA :
F.numFaces = 4;
F.dimension = 2;
- F.monomials = F.serendip ? generatePascalQuadSerendip(F.order) : generatePascalQuad(F.order);
+ F.monomials = F.serendip ? generatePascalQuadSerendip(F.order) :
+ generatePascalQuad(F.order);
F.points = gmshGeneratePointsQuad(F.order, F.serendip);
if (F.order == 0) {
generateClosureOrder0(F.closures, 8);
generateClosureOrder0(F.fullClosures, 8);
F.closureRef.resize(8, 0);
- } else {
+ }
+ else {
generate2dEdgeClosure(F.closures, F.order, 4);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, F.order, 4, F.serendip);
}
@@ -1591,13 +1573,15 @@ const polynomialBasis *polynomialBases::find(int tag)
case TYPE_TET :
F.numFaces = 4;
F.dimension = 3;
- F.monomials = F.serendip ? generatePascalSerendipityTetrahedron(F.order) : generatePascalTetrahedron(F.order);
+ F.monomials = F.serendip ? generatePascalSerendipityTetrahedron(F.order) :
+ generatePascalTetrahedron(F.order);
F.points = gmshGeneratePointsTetrahedron(F.order, F.serendip);
if (F.order == 0) {
generateClosureOrder0(F.closures,24);
generateClosureOrder0(F.fullClosures, 24);
F.closureRef.resize(24, 0);
- } else {
+ }
+ else {
generateFaceClosureTet(F.closures, F.order);
generateFaceClosureTetFull(F.fullClosures, F.closureRef, F.order, F.serendip);
}
@@ -1611,7 +1595,8 @@ const polynomialBasis *polynomialBases::find(int tag)
generateClosureOrder0(F.closures,48);
generateClosureOrder0(F.fullClosures,48);
F.closureRef.resize(48, 0);
- } else {
+ }
+ else {
generateFaceClosurePrism(F.closures, F.order);
generateFaceClosurePrismFull(F.fullClosures, F.closureRef, F.order);
}
@@ -1621,8 +1606,8 @@ const polynomialBasis *polynomialBases::find(int tag)
F.dimension = 3;
F.monomials = generatePascalHex(F.order, F.serendip);
F.points = gmshGeneratePointsHex(F.order, F.serendip);
- // generateFaceClosureHex(F.closures, F.order);
- // generateFaceClosureHexFull(F.fullClosures, F.closureRef, F.order, F.serendip);
+ // generateFaceClosureHex(F.closures, F.order);
+ // generateFaceClosureHexFull(F.fullClosures, F.closureRef, F.order, F.serendip);
break;
}
F.coefficients = generateLagrangeMonomialCoefficients(F.monomials, F.points);
@@ -1655,4 +1640,3 @@ const fullMatrix<double> &polynomialBases::findInjector(int tag1, int tag2)
injector.insert(std::make_pair(key, inj));
return injector[key];
}
-
--
GitLab