diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 9e2390629535053bbac9bcd13fe8fe6e61ab75c5..bccac6128cf1806a20507a7bf9a5fa79317f6865 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -263,8 +263,7 @@ double MElement::getJacobianDeterminant(double u, double v, double w) {
return getJacobian(u,v,w,jac);
}
-
-double MElement::getJacobian(double gsf[][3], double jac[3][3])
+double MElement::getJacobian(const fullMatrix<double> &gsf, double jac[3][3])
{
jac[0][0] = jac[0][1] = jac[0][2] = 0.;
jac[1][0] = jac[1][1] = jac[1][2] = 0.;
@@ -272,11 +271,10 @@ double MElement::getJacobian(double gsf[][3], double jac[3][3])
for (int i = 0; i < getNumVertices(); i++) {
const MVertex* v = getVertex(i);
- double* gg = gsf[i];
for (int j = 0; j < 3; j++) {
- jac[j][0] += v->x() * gg[j];
- jac[j][1] += v->y() * gg[j];
- jac[j][2] += v->z() * gg[j];
+ jac[j][0] += v->x() * gsf(i, j);
+ jac[j][1] += v->y() * gsf(i, j);
+ jac[j][2] += v->z() * gsf(i, j);
}
}
@@ -1039,6 +1037,64 @@ MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
}
}
+/*const fullMatrix<double> &MElement::getShapeFunctionsAtIntegrationPoints (int integrationOrder, int functionSpaceOrder) {
+ static std::map <std::pair<int,int>, fullMatrix<double> > F;
+ const polynomialBasis *fs = getFunctionSpace (functionSpaceOrder);
+ fullMatrix<double> &mat = F [std::make_pair(fs->type, integrationOrder)];
+ if (mat.size1()!=0) return mat;
+ int npts;
+ IntPt *pts;
+ getIntegrationPoints (integrationOrder, &npts, &pts);
+ mat.resize (fs->points.size1(), npts);
+ double f[512];
+ for (int i = 0; i < npts; i++) {
+ fs->f (pts[i].pt[0], pts[i].pt[1], pts[i].pt[2], f);
+ for (int j = 0; j < fs->points.size1(); j++) {
+ mat (i, j) = f[j];
+ }
+ }
+ return mat;
+}*/
+
+const fullMatrix<double> &MElement::getGradShapeFunctionsAtIntegrationPoints (int integrationOrder, int functionSpaceOrder) {
+ static std::map <std::pair<int,int>, fullMatrix<double> > DF;
+ const polynomialBasis *fs = getFunctionSpace (functionSpaceOrder);
+ fullMatrix<double> &mat = DF [std::make_pair(fs->type, integrationOrder)];
+ if (mat.size1()!=0) return mat;
+ int npts;
+ IntPt *pts;
+ getIntegrationPoints (integrationOrder, &npts, &pts);
+ mat.resize (fs->points.size1(), npts*3);
+ double df[512][3];
+ for (int i = 0; i < npts; i++) {
+ fs->df (pts[i].pt[0], pts[i].pt[1], pts[i].pt[2], df);
+ for (int j = 0; j < fs->points.size1(); j++) {
+ mat (j, i*3+0) = df[j][0];
+ mat (j, i*3+1) = df[j][1];
+ mat (j, i*3+2) = df[j][2];
+ }
+ }
+ return mat;
+}
+const fullMatrix<double> &MElement::getGradShapeFunctionsAtNodes (int functionSpaceOrder) {
+ static std::map <int, fullMatrix<double> > DF;
+ const polynomialBasis *fs = getFunctionSpace (functionSpaceOrder);
+ fullMatrix<double> &mat = DF [fs->type];
+ if (mat.size1()!=0) return mat;
+ const fullMatrix<double> &points = fs->points;
+ mat.resize (points.size1(), points.size1()*3);
+ double df[512][3];
+ for (int i = 0; i < points.size1(); i++) {
+ fs->df (points(i,0), points(i,1), points(i,2), df);
+ for (int j = 0; j < points.size1(); j++) {
+ mat (j, i*3+0) = df[j][0];
+ mat (j, i*3+1) = df[j][1];
+ mat (j, i*3+2) = df[j][2];
+ }
+ }
+ return mat;
+}
+
#include "Bindings.h"
void MElement::registerBindings(binding *b)
diff --git a/Geo/MElement.h b/Geo/MElement.h
index ee26dfce960c6121cb0f42ae4e45de6305f058d0..0d5ac4bc74bcbdb69afb5e33b186823fb4fc8f28 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -88,8 +88,7 @@ class MElement
// get the vertex using the VTK ordering
virtual MVertex *getVertexVTK(int num){ return getVertex(num); }
-
- // get the vertex using the Nastran BDF ordering
+// get the vertex using the Nastran BDF ordering
virtual MVertex *getVertexBDF(int num){ return getVertex(num); }
// get the vertex using MED ordering
@@ -223,9 +222,12 @@ class MElement
int order=-1);
virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3],
int order=-1);
+ //const fullMatrix<double> &getShapeFunctionsAtIntegrationPoints (int integrationOrder, int functionSpaceOrder=-1);
+ const fullMatrix<double> &getGradShapeFunctionsAtIntegrationPoints (int integrationOrder, int functionSpaceOrder=-1);
+ const fullMatrix<double> &getGradShapeFunctionsAtNodes (int functionSpaceOrder=-1);
// return the Jacobian of the element evaluated at point (u,v,w) in
// parametric coordinates
- double getJacobian(double gsf[][3], double jac[3][3]);
+ double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]);
double getJacobian(double u, double v, double w, double jac[3][3]);
double getPrimaryJacobian(double u, double v, double w, double jac[3][3]);
//bindings : double[3][3] is not easy to bind, we could use fullMatrix instead
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index 6ae19d0f941320fb773ff3cae557f509519afc14..fb735a6cc261e24c1dfd57e5d7cea837b13b64c6 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -308,44 +308,6 @@ class MTetrahedron10 : public MTetrahedron {
*
*/
-static int reverseTet20[20] = {0,2,1,3, // principal vertices
- 9,8, // E0 switches with E2
- 7,6, // E1 inverts direction
- 5,4, // E2 switches with E0
- 10,11, // E3 pure w edge > remains the same
- 14,15, // E4 uw edge swithes with v/w edge E5
- 12,13, // E5 switches with E4
- 16, // F0 is uv plane, reverts normal
- 18, // F1 is uw plane, switches with F2
- 17, // F2 is vw plane, switches with F1
- 19}; // F3 is uvw plane, reverts normal
-
-static int reverseTet35[35] = {0,2,1,3, // principal vertices
-
- 12,11,10, // E0 switches with E2
- 9,8,7, // E1 inverts direction
- 6,5,4, // E2 switches with E0
- 13,14,15, // E3 pure w edge > remains the same
- 19,20,21, // E4 uw edge swithes with v/w edge E5
- 16,17,18, // E5 switches with E4
- 22,24,23, // F0 is uv plane, reverts normal
- 28,30,29, // F1 is uw plane, switches with F2, orientation is different
- 25,27,26, // F2 is vw plane, switches with F1
- 31,33,32, // F3 is uvw plane, reverts normal
- 34}; // central node remains
-
-static int reverseTet34[34] = {0,2,1,3, // principal vertices
- 12,11,10, // E0 switches with E2
- 9,8,7, // E1 inverts direction
- 6,5,4, // E2 switches with E0
- 13,14,15, // E3 pure w edge > remains the same
- 19,20,21, // E4 uw edge swithes with v/w edge E5
- 16,17,18, // E5 switches with E4
- 22,24,23, // F0 is uv plane, reverts normal
- 28,29,30, // F1 is uw plane, switches with F2
- 25,26,27, // F2 is vw plane, switches with F1
- 31,33,32};// F3 is uvw plane, reverts normal
-
class MTetrahedronN : public MTetrahedron {
const static std::vector<int> &_getReverseIndices(int order);
protected:
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index d10ff18ba3ccc980486a131d8635ef82745fce74..c95bf5618af3d4af2c5e23ef7b7070fc170e6406 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -27,36 +27,31 @@
static bool mappingIsInvertible(MTetrahedron *e)
{
- if (e->getPolynomialOrder() == 1) return 1.0;
+ if (e->getPolynomialOrder() == 1) return true;
double mat[3][3];
e->getPrimaryJacobian(0., 0., 0., mat);
double det0 = det3x3(mat);
-
- IntPt *pts;
- int npts;
- e->getIntegrationPoints(e->getPolynomialOrder(), &npts, &pts);
- for (int i = 0; i < npts; i++){
- const double u = pts[i].pt[0];
- const double v = pts[i].pt[1];
- const double w = pts[i].pt[2];
- e->getJacobian(u, v, w, mat);
- double detN = det3x3(mat);
- if (det0 * detN <= 0.) return false;
+ fullMatrix<double> df;
+ {
+ const fullMatrix<double> &alldf = e->getGradShapeFunctionsAtIntegrationPoints(e->getPolynomialOrder());
+ for (int i = 0; i < alldf.size2()/3; i++){
+ df.setAsProxy(alldf, 3*i, 3);
+ e->getJacobian(df, mat);
+ if (det0 * det3x3(mat) <= 0.) return false;
+ }
}
-
- const fullMatrix<double> &points = e->getFunctionSpace()->points;
-
- for (int i = 0; i < e->getNumPrimaryVertices(); i++) {
- const double u = points(i,0);
- const double v = points(i,1);
- const double w = points(i,2);
- e->getJacobian(u, v, w, mat);
- double detN = det3x3(mat);
- if (det0 * detN <= 0.) return false;
+ {
+ const fullMatrix<double> &points = e->getFunctionSpace()->points;
+ const fullMatrix<double> &alldf = e->getGradShapeFunctionsAtNodes(e->getPolynomialOrder());
+ double gradShapeFunctions[300][3];
+ for (int i = 0; i < alldf.size2()/3; i++){
+ df.setAsProxy(alldf, 3*i, 3);
+ e->getJacobian(df, mat);
+ if (det0 * det3x3(mat) <= 0.) return false;
+ }
}
-
return true;
}
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 629060dedd1fb5ab1f6438974d3fc755aee3f061..7a6ea597a99faf3ac191465bc7ec59270f012eb9 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -1354,6 +1354,7 @@ const polynomialBasis *polynomialBases::find(int tag)
generate3dFaceClosure(F.closures, 1);
break;
}
+ F.type = tag;
F.coefficients = generateLagrangeMonomialCoefficients(F.monomials, F.points);
// printf("Case: %d coeffs:\n",tag);
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index 805df454d6da187b55d144275ab9aefcfe89d101..3779c6716aa73fcba9622598b5ec0401a4acf0d1 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -70,6 +70,8 @@ class binding;
class polynomialBasis
{
public:
+ //for now the only implemented polynomial basis are nodal poly basis, we use the type of the corresponding gmsh element as type
+ int type;
typedef std::vector<std::vector<int> > clCont;
clCont closures;
fullMatrix<double> points;