From 8d0d46e0ca7de6f6ce589afe5e92c83e6952d427 Mon Sep 17 00:00:00 2001
From: Amaury Johnan <amjohnen@gmail.com>
Date: Fri, 1 Aug 2014 15:15:21 +0000
Subject: [PATCH] fix bugs pyramids
---
Numeric/JacobianBasis.cpp | 8 ++++----
Numeric/bezierBasis.cpp | 9 ++++++---
Numeric/pointsGenerators.cpp | 12 +++++++-----
Numeric/pointsGenerators.h | 10 ++++------
4 files changed, 21 insertions(+), 18 deletions(-)
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index 0d1b868377..05ac767ff3 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -54,7 +54,7 @@ GradientBasis::GradientBasis(int tag, int order)
samplingPoints = gmshGeneratePointsHexahedron(order,false);
break;
case TYPE_PYR :
- samplingPoints = JacobianBasis::generateJacPointsPyramid(order);
+ samplingPoints = gmshGeneratePointsPyramidGeneral(false, order+2, order);
break;
default :
Msg::Error("Unknown Jacobian function space for element tag %d", tag);
@@ -92,10 +92,10 @@ void GradientBasis::getGradientsFromNodes(const fullMatrix<double> &nodes,
JacobianBasis::JacobianBasis(int tag, int jacOrder) :
_bezier(NULL), _tag(tag), _dim(ElementType::DimensionFromTag(tag)),
- _jacOrder(jacOrder >= 0 ? jacOrder : JacobianBasis::jacobianOrder(tag))
+ _jacOrder(jacOrder >= 0 ? jacOrder : jacobianOrder(tag))
{
const int parentType = ElementType::ParentTypeFromTag(tag);
- const int primJacobianOrder = JacobianBasis::jacobianOrder(parentType, 1);
+ const int primJacobianOrder = jacobianOrder(parentType, 1);
fullMatrix<double> lagPoints; // Sampling points
@@ -122,7 +122,7 @@ JacobianBasis::JacobianBasis(int tag, int jacOrder) :
lagPoints = gmshGeneratePointsHexahedron(_jacOrder,false);
break;
case TYPE_PYR :
- lagPoints = generateJacPointsPyramid(_jacOrder);
+ lagPoints = gmshGeneratePointsPyramidGeneral(false, _jacOrder+2, _jacOrder);
break;
default :
Msg::Error("Unknown Jacobian function space for element tag %d", tag);
diff --git a/Numeric/bezierBasis.cpp b/Numeric/bezierBasis.cpp
index 92a25a4116..fc8a48f060 100644
--- a/Numeric/bezierBasis.cpp
+++ b/Numeric/bezierBasis.cpp
@@ -252,6 +252,9 @@ std::vector< fullMatrix<double> > generateSubPointsPyr(int order)
fullMatrix<double> ref, prox;
subPoints[0] = gmshGenerateMonomialsPyramidGeneral(false, order+2, order);
+ prox.setAsProxy(subPoints[0], 2, 1);
+ prox.scale(-1);
+ prox.add(order);
subPoints[0].scale(.5/(order+2));
subPoints[1].copy(subPoints[0]);
@@ -380,10 +383,10 @@ fullMatrix<double> generateBez2LagMatrixPyramid
* nChoosek(order, exponent(j, 2))
* pow_int(point(i, 0) / denom, exponent(j, 0))
* pow_int(point(i, 1) / denom, exponent(j, 1))
- * pow_int(point(i, 2) , exponent(j, 2))
+ * pow_int(1. - point(i, 2) , exponent(j, 2))
* pow_int(1. - point(i, 0) / denom, order + 2 - exponent(j, 0))
* pow_int(1. - point(i, 1) / denom, order + 2 - exponent(j, 1))
- * pow_int(1. - point(i, 2) , order - exponent(j, 2));
+ * pow_int(point(i, 2) , order - exponent(j, 2));
}
}
return bez2Lag;
@@ -527,7 +530,7 @@ void bezierBasis::_construct(int parentType, int p)
bezierPoints.resize(_exponents.size1(), _exponents.size2());
const double ord = _order + 2;
for (int i = 0; i < bezierPoints.size1(); ++i) {
- bezierPoints(i, 2) = _exponents(i, 2) / ord;
+ bezierPoints(i, 2) = (_order - _exponents(i, 2)) / ord;
const double scale = 1. - bezierPoints(i, 2);
bezierPoints(i, 0) = _exponents(i, 0) / ord * scale;
bezierPoints(i, 1) = _exponents(i, 1) / ord * scale;
diff --git a/Numeric/pointsGenerators.cpp b/Numeric/pointsGenerators.cpp
index 485bf98575..69ac1fed18 100644
--- a/Numeric/pointsGenerators.cpp
+++ b/Numeric/pointsGenerators.cpp
@@ -107,9 +107,9 @@ fullMatrix<double> gmshGeneratePointsPyramidGeneral(bool pyr, int nij, int nk, b
for (int i = 0; i < points.size1(); ++i) {
int div = pyr ? nk+nij : std::max(nij, nk);
points(i, 2) = (nk - points(i, 2)) / div;
- const double duv = -1. + points(i, 2);
- points(i, 0) = duv + points(i, 0) * 2. / div;
- points(i, 1) = duv + points(i, 1) * 2. / div;
+ const double scale = 1. - points(i, 2);
+ points(i, 0) = scale * (-1 + points(i, 0) * 2. / div);
+ points(i, 1) = scale * (-1 + points(i, 1) * 2. / div);
}
return points;
}
@@ -833,6 +833,8 @@ fullMatrix<double> gmshGenerateMonomialsPyramidGeneral(
monomials(3, 1) = nijBase;
monomials(3, 2) = nk;
+ int index = 4;
+
if (nk > 0) {
monomials(4, 0) = 0;
monomials(4, 1) = 0;
@@ -849,9 +851,9 @@ fullMatrix<double> gmshGenerateMonomialsPyramidGeneral(
monomials(7, 0) = 0;
monomials(7, 1) = nij;
monomials(7, 2) = 0;
- }
- int index = 8;
+ index = 8;
+ }
// Base
if (nijBase > 1) {
diff --git a/Numeric/pointsGenerators.h b/Numeric/pointsGenerators.h
index c3f88f39ab..2d581f2f9d 100644
--- a/Numeric/pointsGenerators.h
+++ b/Numeric/pointsGenerators.h
@@ -28,10 +28,6 @@ fullMatrix<double> gmshGeneratePointsHexahedron(int order, bool serendip = false
fullMatrix<double> gmshGeneratePointsPrism(int order, bool serendip = false);
fullMatrix<double> gmshGeneratePointsPyramid(int order, bool serendip = false);
-
-// If 'serendip' == true, generate points for serendipity pyramid of order 'nk',
-// else if 'pyr' == true, generate points for space {X^i Y^j Z^k | i,j <= k+'nij', k <= 'nk'},
-// else if 'pyr' == false, generate points for space {X^i Y^j Z^k | i,j <= 'nij', k <= 'nk'}
fullMatrix<double> gmshGeneratePointsPyramidGeneral(bool pyr, int nij, int nk, bool serendip = false);
@@ -49,9 +45,11 @@ fullMatrix<double> gmshGenerateMonomialsHexaSerendipity(int order);
fullMatrix<double> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints = false);
fullMatrix<double> gmshGenerateMonomialsPrismSerendipity(int order);
+// Generate monomials of pyramidal nodal space {X^i Y^j Z^k | i,j <= k, k <= 'order'},
fullMatrix<double> gmshGenerateMonomialsPyramid(int order, bool forSerendipPoints = false);
-
-// See explanations for 'gmshGeneratePointsPyramidGeneral'
+// If 'serendip' == true, generate monomials of serendipity pyramid at order 'nk',
+// else if 'pyr' == true, generate monomials of space {X^i Y^j Z^k | i,j <= k+'nij', k <= 'nk'},
+// else if 'pyr' == false, generate monomials of space {X^i Y^j Z^k | i,j <= 'nij', k <= 'nk'}
fullMatrix<double> gmshGenerateMonomialsPyramidGeneral(bool pyr, int nij, int nk, bool forSerendipPoints = false);
--
GitLab