From 2b6f6ac02e3e6ea100b3b0a0dd3dc0a76d043107 Mon Sep 17 00:00:00 2001
From: Amaury Johnan <amjohnen@gmail.com>
Date: Mon, 17 Mar 2014 10:41:08 +0000
Subject: [PATCH]
---
Geo/MPyramid.h | 1 -
Numeric/JacobianBasis.cpp | 18 +++++++++++++++---
Numeric/JacobianBasis.h | 10 ++++++++--
Numeric/bezierBasis.cpp | 34 +++++++++++++++++++++++-----------
4 files changed, 46 insertions(+), 17 deletions(-)
diff --git a/Geo/MPyramid.h b/Geo/MPyramid.h
index 7a3623ffab..7bfdc27caf 100644
--- a/Geo/MPyramid.h
+++ b/Geo/MPyramid.h
@@ -165,7 +165,6 @@ class MPyramid : public MElement {
}
static int faces_pyramid(const int face, const int vert)
{
- // only triangular faces
static const int f[5][4] = {
{0, 1, 4, -1},
{3, 0, 4, -1},
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index af719e909c..20af53c78f 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -27,8 +27,8 @@ JacobianBasis::JacobianBasis(int tag)
{
const int parentType = ElementType::ParentTypeFromTag(tag);
const int order = ElementType::OrderFromTag(tag);
- const int jacobianOrder = JacobianBasis::jacobianOrder(parentType, order);
- const int primJacobianOrder = JacobianBasis::jacobianOrder(parentType, 1);
+ const int jacobianOrder = JacobianBasis::jacobianOrder(parentType, order);
+ const int primJacobianOrder = JacobianBasis::jacobianOrder(parentType, 1);
fullMatrix<double> lagPoints; // Sampling points
@@ -728,7 +728,8 @@ fullMatrix<double> JacobianBasis::generateJacPointsPyramid(int order)
return points;
}
-int JacobianBasis::jacobianOrder(int parentType, int order) {
+int JacobianBasis::jacobianOrder(int parentType, int order)
+{
switch (parentType) {
case TYPE_PNT : return 0;
case TYPE_LIN : return order - 1;
@@ -745,3 +746,14 @@ int JacobianBasis::jacobianOrder(int parentType, int order) {
return 0;
}
}
+
+
+void JacobianBasis::getGradientsFromNodes(const fullMatrix<double> &nodes,
+ fullMatrix<double> *dxyzdX,
+ fullMatrix<double> *dxyzdY,
+ fullMatrix<double> *dxyzdZ) const
+{
+ if (dxyzdX) gradShapeMatX.mult(nodes, *dxyzdX);
+ if (dxyzdY) gradShapeMatY.mult(nodes, *dxyzdY);
+ if (dxyzdZ) gradShapeMatZ.mult(nodes, *dxyzdZ);
+}
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index f1655c2e16..8a08f18e69 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -14,8 +14,6 @@
class JacobianBasis {
private:
- const bezierBasis *bezier;
-
fullMatrix<double> gradShapeMatX, gradShapeMatY, gradShapeMatZ;
fullMatrix<double> gradShapeMatXFast, gradShapeMatYFast, gradShapeMatZFast;
fullVector<double> primGradShapeBarycenterX, primGradShapeBarycenterY, primGradShapeBarycenterZ;
@@ -41,6 +39,8 @@ class JacobianBasis {
fullMatrix<double> &JDJ) const;
public :
+ const bezierBasis *bezier;
+
JacobianBasis(int tag);
// Get methods
@@ -106,6 +106,12 @@ class JacobianBasis {
}
// Jacobian basis order and pyramidal basis
+ void getGradientsFromNodes(const fullMatrix<double> &nodes,
+ fullMatrix<double> *dxyzdX,
+ fullMatrix<double> *dxyzdY,
+ fullMatrix<double> *dxyzdZ) const;
+
+ //
static int jacobianOrder(int parentType, int order);
static fullMatrix<double> generateJacMonomialsPyramid(int order);
static fullMatrix<double> generateJacPointsPyramid(int order);
diff --git a/Numeric/bezierBasis.cpp b/Numeric/bezierBasis.cpp
index ae93f04588..054fdaa958 100644
--- a/Numeric/bezierBasis.cpp
+++ b/Numeric/bezierBasis.cpp
@@ -247,6 +247,10 @@ namespace {
return subPoints;
}
+ Msg::Error("Subpoints for pyramid p>1 not implemented !");
+ return std::vector< fullMatrix<double> >();
+
+ /*
std::vector< fullMatrix<double> > subPoints(8);
fullMatrix<double> ref, prox;
@@ -296,6 +300,7 @@ namespace {
}
return subPoints;
+ */
}
// Matrices generation
@@ -378,18 +383,18 @@ namespace {
for (int j = 0; j < ndofs; j++) {
bez2Lag(i, j) =
nChoosek(order, exponent(j, 2))
- * pow(point(i, 2), exponent(j, 2)) //FIXME use pow_int
- * pow(1. - point(i, 2), order - exponent(j, 2));
-
+ * pow_int(point(i, 2), exponent(j, 2))
+ * pow_int(1. - point(i, 2), order - exponent(j, 2));
+
double p = order + 2 - exponent(j, 2);
double denom = 1. - point(i, 2);
bez2Lag(i, j) *=
nChoosek(p, exponent(j, 0))
* nChoosek(p, exponent(j, 1))
- * pow(point(i, 0) / denom, exponent(j, 0))
- * pow(point(i, 1) / denom, exponent(j, 1))
- * pow(1. - point(i, 0) / denom, p - exponent(j, 0))
- * pow(1. - point(i, 1) / denom, p - exponent(j, 1));
+ * pow_int(point(i, 0) / denom, exponent(j, 0))
+ * pow_int(point(i, 1) / denom, exponent(j, 1))
+ * pow_int(1. - point(i, 0) / denom, p - exponent(j, 0))
+ * pow_int(1. - point(i, 1) / denom, p - exponent(j, 1));
}
}
return bez2Lag;
@@ -458,16 +463,23 @@ void bezierBasis::_construct(int parentType, int p)
dim = 3;
numLagCoeff = 8;
exponents = JacobianBasis::generateJacMonomialsPyramid(order);
- subPoints = generateSubPointsPyr(order);
-
- numDivisions = static_cast<int>(subPoints.size());
+ if (order < 2) {
+ subPoints = generateSubPointsPyr(order);
+ numDivisions = static_cast<int>(subPoints.size());
+ }
+ else {
+ numDivisions = 8;
+ static int a = 0;
+ if (++a == 1) Msg::Warning("Subdivision not available for pyramid p>1");
+ }
fullMatrix<double> bezierPoints = exponents;
bezierPoints.scale(1. / (order + 2));
matrixBez2Lag = generateBez2LagMatrixPyramid(exponents, bezierPoints, order);
matrixBez2Lag.invert(matrixLag2Bez);
- subDivisor = generateSubDivisorPyramid(exponents, subPoints, matrixLag2Bez, order);
+ if (order < 2)
+ subDivisor = generateSubDivisorPyramid(exponents, subPoints, matrixLag2Bez, order);
return;
}
--
GitLab