Skip to content
Snippets Groups Projects
Commit 061311ce authored by Amaury Johnen's avatar Amaury Johnen
Browse files

new algo for pyramids. Validation points ok

parent a7596338
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Common/Gmsh.cpp
+ 5
− 3
View file @ 061311ce
......@@ -244,7 +244,7 @@ int GmshBatch()
solver_batch_cb(0, (void*)CTX::instance()->launchSolverAtStartup);
#endif
if (false) {
if (true) {
// 11/06/13 : This will be removed later !
for (int i = 1; i < MSH_NUM_TYPE+1; ++i) {
if (i == 76 || i == 77 || i == 78)
......@@ -258,9 +258,11 @@ int GmshBatch()
continue;
}
if (MElement::ParentTypeFromTag(i) == TYPE_PYR) {
Msg::Warning("ignoring pyramid %d", i);
if (i > 124 && i < 132) {
Msg::Warning("ignoring serendipity pyramid %d (bug in Bergot implementation)", i);
continue;
}
}
if (MElement::ParentTypeFromTag(i) == TYPE_POLYG) {
Msg::Warning("ignoring polygone %d", i);
continue;
......
Geo/MElement.cpp
+ 14
− 7
View file @ 061311ce
......@@ -1282,11 +1282,18 @@ int MElement::getInfoMSH(const int typeMSH, const char **const name)
case MSH_PYR_14 : if(name) *name = "Pyramid 14"; return 5 + 8 + 1;
case MSH_PYR_30 : if(name) *name = "Pyramid 30"; return 5 + 8*2 + 4*1 + 1*4 + 1;
case MSH_PYR_55 : if(name) *name = "Pyramid 55"; return 5 + 8*3 + 4*3 + 1*9 + 5;
case MSH_PYR_91 : if(name) *name = "Pyramid 91"; return 91;
case MSH_PYR_140 : if(name) *name = "Pyramid 140"; return 140;
case MSH_PYR_204 : if(name) *name = "Pyramid 204"; return 204;
case MSH_PYR_285 : if(name) *name = "Pyramid 285"; return 285;
case MSH_PYR_385 : if(name) *name = "Pyramid 385"; return 385;
case MSH_PYR_91 : if(name) *name = "Pyramid 91"; return 5 + 8*4 + 4*6 + 1*16 + 14;
case MSH_PYR_140 : if(name) *name = "Pyramid 140"; return 5 + 8*5 + 4*10 + 1*25 + 30;
case MSH_PYR_204 : if(name) *name = "Pyramid 204"; return 5 + 8*6 + 4*15 + 1*36 + 55;
case MSH_PYR_285 : if(name) *name = "Pyramid 285"; return 5 + 8*7 + 4*21 + 1*49 + 91;
case MSH_PYR_385 : if(name) *name = "Pyramid 385"; return 5 + 8*8 + 4*28 + 1*64 + 140;
case MSH_PYR_29 : if(name) *name = "Pyramid 29"; return 5 + 8*2 + 4*1 + 1*4;
case MSH_PYR_50 : if(name) *name = "Pyramid 50"; return 5 + 8*3 + 4*3 + 1*9;
case MSH_PYR_77 : if(name) *name = "Pyramid 77"; return 5 + 8*4 + 4*6 + 1*16;
case MSH_PYR_110 : if(name) *name = "Pyramid 110"; return 5 + 8*5 + 4*10 + 1*25;
case MSH_PYR_149 : if(name) *name = "Pyramid 149"; return 5 + 8*6 + 4*15 + 1*36;
case MSH_PYR_194 : if(name) *name = "Pyramid 194"; return 5 + 8*7 + 4*21 + 1*49;
case MSH_PYR_245 : if(name) *name = "Pyramid 245"; return 5 + 8*8 + 4*28 + 1*64;
case MSH_POLYH_ : if(name) *name = "Polyhedron"; return 0;
case MSH_PNT_SUB : if(name) *name = "Point Xfem"; return 1;
case MSH_LIN_SUB : if(name) *name = "Line Xfem"; return 2;
......
Geo/MPyramid.h
+ 6
− 5
View file @ 061311ce
......@@ -166,11 +166,12 @@ class MPyramid : public MElement {
static int faces_pyramid(const int face, const int vert)
{
// only triangular faces
static const int f[4][3] = {
{0, 1, 4},
{3, 0, 4},
{1, 2, 4},
{2, 3, 4}
static const int f[5][4] = {
{0, 1, 4, -1},
{3, 0, 4, -1},
{1, 2, 4, -1},
{2, 3, 4, -1},
{0, 3, 2, 1}
};
return f[face][vert];
}
......
Numeric/nodalBasis.cpp
+ 1
− 1
View file @ 061311ce
......@@ -1606,7 +1606,7 @@ nodalBasis::nodalBasis(int tag)
closureRef.resize(48, 0);
}
else {
if (order < 2) generateFaceClosurePrism(closures, order);
generateFaceClosurePrism(closures, order);
generateFaceClosurePrismFull(fullClosures, closureRef, order);
}
break;
......
Numeric/pointsGenerators.cpp
+ 103
− 8
View file @ 061311ce
......@@ -9,6 +9,7 @@
#include "MTetrahedron.h"
#include "MHexahedron.h"
#include "MPrism.h"
#include "MPyramid.h"
/* --- Lines --- */
......@@ -962,7 +963,8 @@ fullMatrix<int> gmshGenerateMonomialsTetrahedron(int order, bool serendip)
fullMatrix<int> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints)
{
int nbMonomials = forSerendipPoints ? 6 + (order-1)*9 : (order + 1) * (order + 1)*(order + 2)/2;
int nbMonomials = forSerendipPoints ? 6 + (order-1)*9 :
(order + 1) * (order + 1)*(order + 2)/2;
fullMatrix<int> monomials(nbMonomials, 3);
monomials(0, 0) = 0;
......@@ -1028,9 +1030,7 @@ fullMatrix<int> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints)
i2 = MPrism::faces_prism(iface, 2);
dudv.setAsProxy(dudvT);
}
else {
continue;
}
else continue;
int u[3];
u[0] = (monomials(i1, 0) - monomials(i0, 0)) / order;
......@@ -1068,7 +1068,8 @@ fullMatrix<int> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints)
fullMatrix<int> gmshGenerateMonomialsHexahedron(int order, bool forSerendipPoints)
{
int nbMonomials = forSerendipPoints ? 8 + (order-1)*12 : (order+1)*(order+1)*(order+1);
int nbMonomials = forSerendipPoints ? 8 + (order-1)*12 :
(order+1)*(order+1)*(order+1);
fullMatrix<int> monomials(nbMonomials, 3);
monomials(0, 0) = 0;
......@@ -1121,10 +1122,10 @@ fullMatrix<int> gmshGenerateMonomialsHexahedron(int order, bool forSerendipPoint
}
}
if (!forSerendipPoints) {
fullMatrix<int> dudv = gmshGenerateMonomialsQuadrangle(order - 2);
dudv.add(1);
if (!forSerendipPoints) {
for (int iface = 0; iface < 6; ++iface) {
int i0 = MHexahedron::faces_hexa(iface, 0);
int i1 = MHexahedron::faces_hexa(iface, 1);
......@@ -1155,3 +1156,97 @@ fullMatrix<int> gmshGenerateMonomialsHexahedron(int order, bool forSerendipPoint
return monomials;
}
fullMatrix<int> gmshGenerateMonomialsPyramid(int order, bool forSerendipPoints)
{
int nbMonomials = forSerendipPoints ? 5 + (order-1)*8 :
(order+1)*((order+1)+1)*(2*(order+1)+1)/6;
fullMatrix<int> monomials(nbMonomials, 3);
monomials(0, 0) = 0;
monomials(0, 1) = 0;
monomials(0, 2) = 0;
if (order > 0) {
monomials(1, 0) = order;
monomials(1, 1) = 0;
monomials(1, 2) = 0;
monomials(2, 0) = order;
monomials(2, 1) = order;
monomials(2, 2) = 0;
monomials(3, 0) = 0;
monomials(3, 1) = order;
monomials(3, 2) = 0;
monomials(4, 0) = 0;
monomials(4, 1) = 0;
monomials(4, 2) = order;
if (order > 1) {
int index = 5;
for (int iedge = 0; iedge < 8; ++iedge) {
int i0 = MPyramid::edges_pyramid(iedge, 0);
int i1 = MPyramid::edges_pyramid(iedge, 1);
int u_1 = (monomials(i1,0)-monomials(i0,0)) / order;
int u_2 = (monomials(i1,1)-monomials(i0,1)) / order;
int u_3 = (monomials(i1,2)-monomials(i0,2)) / order;
for (int i = 1; i < order; ++i, ++index) {
monomials(index, 0) = monomials(i0, 0) + i * u_1;
monomials(index, 1) = monomials(i0, 1) + i * u_2;
monomials(index, 2) = monomials(i0, 2) + i * u_3;
}
}
if (!forSerendipPoints) {
fullMatrix<int> dudvQ = gmshGenerateMonomialsQuadrangle(order - 2);
dudvQ.add(1);
fullMatrix<int> dudvT;
if (order > 2) dudvT = gmshGenerateMonomialsTriangle(order - 3);
dudvT.add(1);
for (int iface = 0; iface < 5; ++iface) {
int i0, i1, i2;
i0 = MPyramid::faces_pyramid(iface, 0);
i1 = MPyramid::faces_pyramid(iface, 1);
fullMatrix<int> dudv;
if (MPyramid::faces_pyramid(iface, 3) != -1) {
i2 = MPyramid::faces_pyramid(iface, 3);
dudv.setAsProxy(dudvQ);
}
else if (order > 2) {
i2 = MPyramid::faces_pyramid(iface, 2);
dudv.setAsProxy(dudvT);
}
else continue;
int u[3];
u[0] = (monomials(i1, 0) - monomials(i0, 0)) / order;
u[1] = (monomials(i1, 1) - monomials(i0, 1)) / order;
u[2] = (monomials(i1, 2) - monomials(i0, 2)) / order;
int v[3];
v[0] = (monomials(i2, 0) - monomials(i0, 0)) / order;
v[1] = (monomials(i2, 1) - monomials(i0, 1)) / order;
v[2] = (monomials(i2, 2) - monomials(i0, 2)) / order;
for (int i = 0; i < dudv.size1(); ++i, ++index) {
monomials(index, 0) = monomials(i0, 0) + u[0] * dudv(i, 0) + v[0] * dudv(i, 1);
monomials(index, 1) = monomials(i0, 1) + u[1] * dudv(i, 0) + v[1] * dudv(i, 1);
monomials(index, 2) = monomials(i0, 2) + u[2] * dudv(i, 0) + v[2] * dudv(i, 1);
}
}
if (order > 2) {
fullMatrix<int> inner = gmshGenerateMonomialsPyramid(order - 3);
inner.add(1);
monomials.copy(inner, 0, nbMonomials - index, 0, 3, index, 0);
}
}
}
}
return monomials;
}
Numeric/pointsGenerators.h
+ 2
− 1
View file @ 061311ce
......@@ -44,7 +44,8 @@ fullMatrix<int> gmshGenerateMonomialsHexahedron(int order, bool forSerendipPoint
fullMatrix<int> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints = false);
//fullMatrix<int> gmshGenerateMonomialsPrismSerendipity(int order); //FIXME
//fullMatrix<int> gmshGenerateMonomialsPyramid(int order, bool serendip);
fullMatrix<int> gmshGenerateMonomialsPyramid(int order, bool forSerendipPoints = false);
//fullMatrix<int> gmshGenerateMonomialsPyramidSerendipity(int order); //FIXME
......
Numeric/polynomialBasis.h
+ 0
− 2
View file @ 061311ce
......@@ -75,8 +75,6 @@ class polynomialBasis : public nodalBasis
polynomialBasis(int tag);
~polynomialBasis();
//int compareNewAlgoPointsWithOld() const;
virtual inline int getNumShapeFunctions() const {return coefficients.size1();}
virtual void f(double u, double v, double w, double *sf) const;
......
Numeric/pyramidalBasis.cpp
+ 26
− 2
View file @ 061311ce
......@@ -7,6 +7,15 @@
#include "pointsGenerators.h"
static void copy(const fullMatrix<int> &orig, fullMatrix<double> &b)
{
b.resize(orig.size1(), orig.size2());
for (int i = 0; i < orig.size1(); ++i) {
for (int j = 0; j < orig.size2(); ++j) {
b(i, j) = static_cast<double>(orig(i, j));
}
}
}
pyramidalBasis::pyramidalBasis(int tag) : nodalBasis(tag)
{
......@@ -28,6 +37,21 @@ pyramidalBasis::pyramidalBasis(int tag) : nodalBasis(tag)
delete[] fval;
// TEST NEW ALGO POINTS / MONOMIAL
monomials_newAlgo = gmshGenerateMonomialsPyramid(order, serendip);
copy(monomials_newAlgo, points_newAlgo);
if (order == 0) return;
for (int i = 0; i < points_newAlgo.size1(); ++i) {
points_newAlgo(i, 2) = points_newAlgo(i, 2) / order;
if (i != 4) {
const double duv = -1. + points_newAlgo(i, 2);
points_newAlgo(i, 0) = duv + points_newAlgo(i, 0) * 2. / order;
points_newAlgo(i, 1) = duv + points_newAlgo(i, 1) * 2. / order;
}
}
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment