diff --git a/CMakeLists.txt b/CMakeLists.txt
index e55427ad08b5a84b22483be83e41ec543fbcfb13..628a5086223ee51298a5733a1eea18cb252d101c 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -672,13 +672,16 @@ if(ENABLE_OCC)
${OCC_SYS_NAME}/lib)
if(OCC_LIB)
list(APPEND OCC_LIBS ${OCC_LIB})
+ else(OCC_LIB)
+ message(STATUS "OCC lib " ${OCC} " not Found")
endif(OCC_LIB)
set(OCC_LIB OCC_LIB-NOTFOUND CACHE INTERNAL "")
# unset(OCC_LIB CACHE) # cleaner, but only available in cmake >= 2.6.4
endforeach(OCC)
list(LENGTH OCC_LIBS NUM_OCC_LIBS)
- if(NUM_OCC_LIBS EQUAL NUM_OCC_LIBS_REQUIRED)
+
+ if(NUM_OCC_LIBS EQUAL NUM_OCC_LIBS_REQUIRED)
find_path(OCC_INC "BRep_Tool.hxx" PATHS ENV CASROOT PATH_SUFFIXES inc
include opencascade)
if(OCC_INC)
diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index 3c23232d10bef3c2e1ef2f8067b9925f5c7d1941..60ad590ff09f66c0b5d8a78e1f584f86c5a4e802 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -134,9 +134,6 @@
#define MSH_TET_165 73
#define MSH_TET_220 74
#define MSH_TET_286 75
-#define MSH_HEX_64 76
-#define MSH_HEX_125 77
-#define MSH_HEX_196 78
#define MSH_TET_74 79
#define MSH_TET_100 80
#define MSH_TET_130 81
@@ -150,7 +147,22 @@
#define MSH_PRI_1 89
#define MSH_PRI_40 90
#define MSH_PRI_75 91
-#define MSH_PRI_126 92
+// HEXES COMPLETE (3->9)
+#define MSH_HEX_64 92
+#define MSH_HEX_125 93
+#define MSH_HEX_216 94
+#define MSH_HEX_343 95
+#define MSH_HEX_512 96
+#define MSH_HEX_729 97
+#define MSH_HEX_1000 98
+// HEXES INCOMPLETE (3->9)
+#define MSH_HEX_56 99
+#define MSH_HEX_98 100
+#define MSH_HEX_152 101
+#define MSH_HEX_222 102
+#define MSH_HEX_296 103
+#define MSH_HEX_386 104
+#define MSH_HEX_488 105
#define MSH_NUM_TYPE 92
diff --git a/Common/GmshMessage.cpp b/Common/GmshMessage.cpp
index a7c6357060fc7cc4dff2fec3eb4b1339ac13f4fa..73faaa148560803ea2149c1d7384aecb72f11923 100644
--- a/Common/GmshMessage.cpp
+++ b/Common/GmshMessage.cpp
@@ -73,6 +73,7 @@ void Msg::Init(int argc, char **argv)
#endif
#if defined(HAVE_PETSC)
PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
+ printf("coucou\n");
#endif
#if defined(HAVE_SLEPC)
SlepcInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 40ca858f718aee2698bb0591896b59441ad5d9a0..5c52edb5b72783283f6a536261e1ba91dfa0730a 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -255,7 +255,7 @@ double MElement::getJacobian(double u, double v, double w, double jac[3][3])
jac[1][0] = jac[1][1] = jac[1][2] = 0.;
jac[2][0] = jac[2][1] = jac[2][2] = 0.;
- double gsf[256][3];
+ double gsf[1256][3];
getGradShapeFunctions(u, v, w, gsf);
for (int i = 0; i < getNumShapeFunctions(); i++) {
const MVertex *v = getShapeFunctionNode(i);
@@ -298,7 +298,7 @@ double MElement::getPrimaryJacobian(double u, double v, double w, double jac[3][
jac[1][0] = jac[1][1] = jac[1][2] = 0.;
jac[2][0] = jac[2][1] = jac[2][2] = 0.;
- double gsf[256][3];
+ double gsf[1256][3];
getGradShapeFunctions(u, v, w, gsf, 1);
for(int i = 0; i < getNumPrimaryShapeFunctions(); i++) {
const MVertex *v = getShapeFunctionNode(i);
@@ -316,7 +316,7 @@ double MElement::getPrimaryJacobian(double u, double v, double w, double jac[3][
void MElement::pnt(double u, double v, double w, SPoint3 &p)
{
double x = 0., y = 0., z = 0.;
- double sf[256];
+ double sf[1256];
getShapeFunctions(u, v, w, sf);
for (int j = 0; j < getNumShapeFunctions(); j++) {
const MVertex *v = getShapeFunctionNode(j);
@@ -330,7 +330,7 @@ void MElement::pnt(double u, double v, double w, SPoint3 &p)
void MElement::primaryPnt(double u, double v, double w, SPoint3 &p)
{
double x = 0., y = 0., z = 0.;
- double sf[256];
+ double sf[1256];
getShapeFunctions(u, v, w, sf, 1);
for (int j = 0; j < getNumPrimaryShapeFunctions(); j++) {
const MVertex *v = getShapeFunctionNode(j);
@@ -354,7 +354,7 @@ void MElement::xyz2uvw(double xyz[3], double uvw[3])
double jac[3][3];
if(!getJacobian(uvw[0], uvw[1], uvw[2], jac)) break;
double xn = 0., yn = 0., zn = 0.;
- double sf[256];
+ double sf[1256];
getShapeFunctions(uvw[0], uvw[1], uvw[2], sf);
for (int i = 0; i < getNumShapeFunctions(); i++) {
MVertex *v = getShapeFunctionNode(i);
@@ -405,7 +405,7 @@ double MElement::interpolate(double val[], double u, double v, double w, int str
{
double sum = 0;
int j = 0;
- double sf[256];
+ double sf[1256];
getShapeFunctions(u, v, w, sf, order);
for(int i = 0; i < getNumShapeFunctions(); i++){
sum += val[j] * sf[i];
@@ -419,7 +419,7 @@ void MElement::interpolateGrad(double val[], double u, double v, double w, doubl
{
double dfdu[3] = {0., 0., 0.};
int j = 0;
- double gsf[256][3];
+ double gsf[1256][3];
getGradShapeFunctions(u, v, w, gsf, order);
for(int i = 0; i < getNumShapeFunctions(); i++){
dfdu[0] += val[j] * gsf[i][0];
@@ -1003,6 +1003,20 @@ int MElement::getInfoMSH(const int typeMSH, const char **const name)
case MSH_HEX_8 : if(name) *name = "Hexahedron 8"; return 8;
case MSH_HEX_20 : if(name) *name = "Hexahedron 20"; return 8 + 12;
case MSH_HEX_27 : if(name) *name = "Hexahedron 27"; return 8 + 12 + 6 + 1;
+ case MSH_HEX_64 : if(name) *name = "Hexahedron 64"; return 64;
+ case MSH_HEX_125 : if(name) *name = "Hexahedron 125"; return 125;
+ case MSH_HEX_216 : if(name) *name = "Hexahedron 216"; return 216;
+ case MSH_HEX_343 : if(name) *name = "Hexahedron 343"; return 343;
+ case MSH_HEX_512 : if(name) *name = "Hexahedron 512"; return 512;
+ case MSH_HEX_729 : if(name) *name = "Hexahedron 729"; return 729;
+ case MSH_HEX_1000 : if(name) *name = "Hexahedron 1000"; return 1000;
+ case MSH_HEX_56 : if(name) *name = "Hexahedron 56"; return 56;
+ case MSH_HEX_98 : if(name) *name = "Hexahedron 98"; return 98;
+ case MSH_HEX_152 : if(name) *name = "Hexahedron 152"; return 152;
+ case MSH_HEX_222 : if(name) *name = "Hexahedron 222"; return 222;
+ case MSH_HEX_296 : if(name) *name = "Hexahedron 296"; return 296;
+ case MSH_HEX_386 : if(name) *name = "Hexahedron 386"; return 386;
+ case MSH_HEX_488 : if(name) *name = "Hexahedron 488"; return 488;
case MSH_PRI_6 : if(name) *name = "Prism 6"; return 6;
case MSH_PRI_15 : if(name) *name = "Prism 15"; return 6 + 9;
case MSH_PRI_18 : if(name) *name = "Prism 18"; return 6 + 9 + 3;
@@ -1164,6 +1178,14 @@ MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
case MSH_TET_220: return new MTetrahedronN(v, 9, num, part);
case MSH_TET_286: return new MTetrahedronN(v, 10, num, part);
case MSH_POLYH_: return new MPolyhedron(v, num, part, owner, parent);
+ case MSH_HEX_56: return new MHexahedronN(v, 3, num, part);
+ case MSH_HEX_64: return new MHexahedronN(v, 3, num, part);
+ case MSH_HEX_125: return new MHexahedronN(v, 4, num, part);
+ case MSH_HEX_216: return new MHexahedronN(v, 5, num, part);
+ case MSH_HEX_343: return new MHexahedronN(v, 6, num, part);
+ case MSH_HEX_512: return new MHexahedronN(v, 7, num, part);
+ case MSH_HEX_729: return new MHexahedronN(v, 8, num, part);
+ case MSH_HEX_1000: return new MHexahedronN(v, 9, num, part);
default: return 0;
}
}
diff --git a/Geo/MHexahedron.cpp b/Geo/MHexahedron.cpp
index 1f64319cf53a8ae6f73eff4748c010b8fc292de9..69c03e2527338abcdbcfc3040ebb79d3d4a211ab 100644
--- a/Geo/MHexahedron.cpp
+++ b/Geo/MHexahedron.cpp
@@ -5,6 +5,8 @@
#include "MHexahedron.h"
#include "Numeric.h"
+#include "Context.h"
+#include "polynomialBasis.h"
int MHexahedron::getVolumeSign()
{
@@ -73,3 +75,87 @@ void MHexahedron::getFaceInfo(const MFace &face, int &ithFace, int &sign, int &r
}
Msg::Error("Could not get face information for hexahedron %d", getNum());
}
+
+void MHexahedronN::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
+{
+ const int numSubEdges = CTX::instance()->mesh.numSubEdges;
+ static double pp[8][3] = {
+ {-1,-1,-1},{1,-1,-1},{1,1,-1},{-1,1,-1},
+ {-1,-1,1},{1,-1,1},{1,1,1},{-1,1,1} };
+ static int ed [12][2] = {
+ {0,1},{0,3},{0,4},{1,2},{1,5},{2,3},
+ {2,6},{3,7},{4,5},{4,7},{5,6},{7,6}
+ };
+ int iEdge = num / numSubEdges;
+ int iSubEdge = num % numSubEdges;
+
+ int iVertex1 = ed [iEdge][0];
+ int iVertex2 = ed [iEdge][1];
+ double t1 = (double) iSubEdge / (double) numSubEdges;
+ double u1 = pp[iVertex1][0] * (1.-t1) + pp[iVertex2][0] * t1;
+ double v1 = pp[iVertex1][1] * (1.-t1) + pp[iVertex2][1] * t1;
+ double w1 = pp[iVertex1][2] * (1.-t1) + pp[iVertex2][2] * t1;
+
+ double t2 = (double) (iSubEdge+1) / (double) numSubEdges;
+ double u2 = pp[iVertex1][0] * (1.-t2) + pp[iVertex2][0] * t2;
+ double v2 = pp[iVertex1][1] * (1.-t2) + pp[iVertex2][1] * t2;
+ double w2 = pp[iVertex1][2] * (1.-t2) + pp[iVertex2][2] * t2;
+
+ SPoint3 pnt1, pnt2;
+ pnt(u1,v1,w1,pnt1);
+ pnt(u2,v2,w2,pnt2);
+ x[0] = pnt1.x(); x[1] = pnt2.x();
+ y[0] = pnt1.y(); y[1] = pnt2.y();
+ z[0] = pnt1.z(); z[1] = pnt2.z();
+
+ // not great, but better than nothing
+ static const int f[6] = {0, 0, 0, 1, 2, 3};
+ n[0] = n[1] = 1 ;
+}
+
+int MHexahedronN::getNumEdgesRep(){
+ return 12 * CTX::instance()->mesh.numSubEdges;
+}
+
+//int MHexaedronN::getNumFacesRep(){
+// return 8 * SQU(CTX::instance()->mesh.numSubEdges);
+//}
+
+const polynomialBasis* MHexahedron::getFunctionSpace(int o) const
+{
+ int order = (o == -1) ? getPolynomialOrder() : o;
+
+ int nv = getNumVolumeVertices();
+
+ if ((nv == 0) && (o == -1)) {
+ switch (order) {
+ case 0: return polynomialBases::find(MSH_HEX_1);
+ case 1: return polynomialBases::find(MSH_HEX_8);
+ case 2: return polynomialBases::find(MSH_HEX_20);
+ case 3: return polynomialBases::find(MSH_HEX_56);
+ case 4: return polynomialBases::find(MSH_HEX_98);
+ case 5: return polynomialBases::find(MSH_HEX_152);
+ case 6: return polynomialBases::find(MSH_HEX_222);
+ case 7: return polynomialBases::find(MSH_HEX_296);
+ case 8: return polynomialBases::find(MSH_HEX_386);
+ case 9: return polynomialBases::find(MSH_HEX_488);
+ default: Msg::Error("Order %d hex function space not implemented", order);
+ }
+ }
+ else {
+ switch (order) {
+ case 0: return polynomialBases::find(MSH_HEX_1);
+ case 1: return polynomialBases::find(MSH_HEX_8);
+ case 2: return polynomialBases::find(MSH_HEX_27);
+ case 3: return polynomialBases::find(MSH_HEX_64);
+ case 4: return polynomialBases::find(MSH_HEX_125);
+ case 5: return polynomialBases::find(MSH_HEX_216);
+ case 6: return polynomialBases::find(MSH_HEX_343);
+ case 7: return polynomialBases::find(MSH_HEX_512);
+ case 8: return polynomialBases::find(MSH_HEX_729);
+ case 9: return polynomialBases::find(MSH_HEX_1000);
+ default: Msg::Error("Order %d hex function space not implemented", order);
+ }
+ }
+ return 0;
+}
diff --git a/Geo/MHexahedron.h b/Geo/MHexahedron.h
index efb616abd9a203ab4efeac1a16d9b28152a09ac3..40a4563bc6b7789536c823c83233c7dcb60e3b7d 100644
--- a/Geo/MHexahedron.h
+++ b/Geo/MHexahedron.h
@@ -57,6 +57,7 @@ class MHexahedron : public MElement {
virtual int getDim() const { return 3; }
virtual int getNumVertices() const { return 8; }
virtual MVertex *getVertex(int num){ return _v[num]; }
+ virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
virtual MVertex *getVertexMED(int num)
{
static const int map[8] = {0, 3, 2, 1, 4, 7, 6, 5};
@@ -403,7 +404,9 @@ class MHexahedron27 : public MHexahedron {
}
virtual int getNumEdgeVertices() const { return 12; }
virtual int getNumFaceVertices() const { return 6; }
- virtual int getNumVolumeVertices() const { return 1; }
+ virtual int getNumVolumeVertices() const {
+ return 8;
+ }
virtual int getNumEdgesRep(){ return 24; }
virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
{
@@ -490,4 +493,141 @@ class MHexahedron27 : public MHexahedron {
}
};
+/*
+ * MHexahedronN
+ *
+ * 3----13----2
+ * |\ |\
+ * |15 24 | 14
+ * 9 \ 20 11 \
+ * | 7----19+---6
+ * |22 | 26 | 23|
+ * 0---+-8----1 |
+ * \ 17 25 \ 18
+ * 10 | 21 12|
+ * \| \|
+ * 4----16----5
+ *
+ * N = 0RDER - 1
+ * 8 := 8 --> 8 + N
+ * 9 := 8 + N + 1 --> 8 + 2 * N
+ * :
+ * 19 := 8 + 11 * N + 1 --> 8 + 12 * N
+ * 20 := 8 + 12 * N + 1 --> 8 + 12 * N + N^2
+ * 21 := 8 + 12 * N + N^2 --> 8 + 12 * N + 2 * N^2
+ * :
+ * 25 := 8 + 12 * N + 5 * N^2 + 1 --> 8 + 12 * N + 6 * N^2
+ * 26 := 8 + 12 * N + 6 * N^2 + 1 --> 8 + 12 * N + 6 * N^2 + N^3
+ *
+ */
+
+
+class MHexahedronN : public MHexahedron {
+ protected:
+ const char _order;
+ std::vector<MVertex*> _vs;
+ public :
+ MHexahedronN(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3,
+ MVertex *v4, MVertex *v5, MVertex *v6, MVertex *v7,
+ const std::vector<MVertex*> &v, char order, int num=0, int part=0)
+ : MHexahedron(v0,v1,v2,v3,v4,v5,v6,v7, num, part), _order(order), _vs(v)
+ {
+ for(unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
+ }
+ MHexahedronN(const std::vector<MVertex*> &v, char order, int num=0, int part=0)
+ : MHexahedron(v[0], v[1], v[2], v[3],v[4], v[5], v[6], v[7], num, part), _order(order)
+ {
+ for(unsigned int i = 8; i < v.size(); i++) _vs.push_back(v[i]);
+ for(unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
+ }
+
+ ~MHexahedronN(){}
+ virtual int getPolynomialOrder() const { return (int)_order; }
+ virtual int getNumVertices() const { return 8 + _vs.size(); }
+ virtual MVertex *getVertex(int num){ return num < 8 ? _v[num] : _vs[num - 8]; }
+ virtual MVertex *getVertexDIFF(int num)
+ {
+ throw;
+ }
+ virtual int getNumEdgeVertices() const { return 12 * (_order - 1); }
+ virtual int getNumFaceVertices() const { return 6 * (_order - 1)*(_order - 1); }
+ virtual int getNumVolumeVertices() const {
+ switch(getTypeForMSH()){
+ case MSH_HEX_27 :
+ case MSH_HEX_64 :
+ case MSH_HEX_125 :
+ case MSH_HEX_216 :
+ case MSH_HEX_343 :
+ case MSH_HEX_512 :
+ case MSH_HEX_729 :
+ case MSH_HEX_1000 :
+ return (_order - 1)*(_order - 1)*(_order - 1);
+ default:
+ return 0;
+ }
+ }
+ virtual int getNumEdgesRep();
+ virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n);
+ virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
+ {
+ v.resize(_order+1);
+ MHexahedron::_getEdgeVertices(num, v);
+ for (int i=0;i<_order-1;i++) v[2+i] = _vs[(_order-1)*num+i];
+ }
+ //virtual int getNumFacesRep();
+ virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
+ {
+ v.resize((_order+1)*(_order+1));
+ MHexahedron::_getFaceVertices(num, v);
+ static const int f[6][4] = {
+ {0, 3, 2, 1},
+ {0, 1, 5, 4},
+ {0, 4, 7, 3},
+ {1, 2, 6, 5},
+ {2, 3, 7, 6},
+ {4, 5, 6, 7}
+ };
+ int count = 4;
+ for (int i=0;i<4;i++){
+ for (int j=0;j<_order-1;j++) v[count++] = _vs[(_order-1)*f[num][i]+j];
+ }
+ for (int i=0;i<(_order+1)*(_order+1);i++){
+ int N = _order - 1;
+ int start = 8 + 12 * N + num * (_order-1)*(_order-1);
+ v[count++] = _vs[start + i];
+ }
+ }
+ virtual int getTypeForMSH() const
+ {
+ // (p+1)^3
+ if(_order == 3 && _vs.size() + 8 == 64 ) return MSH_HEX_64;
+ if(_order == 4 && _vs.size() + 8 == 125) return MSH_HEX_125;
+ if(_order == 5 && _vs.size() + 8 == 216) return MSH_HEX_216;
+ if(_order == 6 && _vs.size() + 8 == 343) return MSH_HEX_343;
+ if(_order == 7 && _vs.size() + 8 == 512) return MSH_HEX_512;
+ if(_order == 8 && _vs.size() + 8 == 729) return MSH_HEX_729;
+ if(_order == 9 && _vs.size() + 8 == 1000) return MSH_HEX_1000;
+ if(_order == 3 && _vs.size() + 8 == 56 ) return MSH_HEX_56;
+ if(_order == 4 && _vs.size() + 8 == 98) return MSH_HEX_98;
+ if(_order == 5 && _vs.size() + 8 == 152) return MSH_HEX_152;
+ if(_order == 6 && _vs.size() + 8 == 222) return MSH_HEX_222;
+ if(_order == 7 && _vs.size() + 8 == 296) return MSH_HEX_296;
+ if(_order == 8 && _vs.size() + 8 == 386) return MSH_HEX_386;
+ if(_order == 9 && _vs.size() + 8 == 488) return MSH_HEX_488;
+ return 0;
+ }
+ virtual void getShapeFunctions(double u, double v, double w, double s[], int o){
+ MElement::getShapeFunctions (u,v,w,s,o);
+ }
+ virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o)
+ {
+ MElement::getGradShapeFunctions (u,v,w,s,o);
+ }
+
+ virtual void revert()
+ {
+ Msg::Error("not done yet reverse hexN");
+ }
+};
+
#endif
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 3882923b207854eda2b973d04cdbb2fc39e9bb4d..026a3f6b719c82ca1d6e8aa554e9ca7126c0f8dc 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -412,7 +412,7 @@ static void getFaceVertices(GFace *gf, MElement *incomplete, MElement *ele,
}
else{
std::vector<MVertex*> &vtcs = faceVertices[face];
- SPoint2 pts[100];
+ SPoint2 pts[1000];
bool reparamOK = true;
if(!linear){
for(int k = 0; k < incomplete->getNumVertices(); k++)
@@ -430,7 +430,7 @@ static void getFaceVertices(GFace *gf, MElement *incomplete, MElement *ele,
}
else{
double X(0), Y(0), Z(0), GUESS[2] = {0, 0};
- double sf[256];
+ double sf[1256];
incomplete->getShapeFunctions(t1, t2, 0, sf);
for (int j = 0; j < incomplete->getNumShapeFunctions(); j++){
MVertex *vt = incomplete->getShapeFunctionNode(j);
@@ -498,6 +498,71 @@ static void reorientTrianglePoints(std::vector<MVertex*> &vtcs, int orientation,
}
}
+static void reorientQuadPoints(std::vector<MVertex*> &vtcs, int orientation,
+ bool swap, int order)
+{
+ int nbPts = vtcs.size();
+ if (nbPts <= 1) return;
+ std::vector<MVertex*> tmp(nbPts);
+
+ // printf("before : ");
+ // for (int i=0;i<vtcs.size();i++)printf("%d ",vtcs[i]->getNum());
+ // printf("\n");
+ int start = 0;
+ while (1){
+ // CORNERS
+ int index = 0;
+ if (order == 0){
+ start++;
+ }
+ else{
+ int i1,i2,i3,i4;
+ if (!swap){
+ if (orientation == 0){i1 = 0; i2=1; i3=2; i4 = 3;}
+ else if (orientation == 1){i1 = 3; i2=0; i3=1; i4 = 2;}
+ else if (orientation == 2){i1 = 2; i2=3; i3=0; i4 = 1;}
+ else if (orientation == 3){i1 = 1; i2=2; i3=3; i4 = 0;}
+ }
+ else{
+ if (orientation == 0){i1 = 0; i2=3; i3=2; i4 = 1;}
+ else if (orientation == 3){i1 = 3; i2=2; i3=1; i4 = 0;}
+ else if (orientation == 2){i1 = 2; i2=1; i3=0; i4 = 3;}
+ else if (orientation == 1){i1 = 1; i2=0; i3=3; i4 = 2;}
+ }
+
+ int indices[4] = {i1,i2,i3,i4};
+ for (int i=0;i<4;i++)tmp[i] = vtcs[start+indices[i]];
+ for (int i=0;i<4;i++)vtcs[start+i] = tmp[i];
+
+ // EDGES
+ for (int iEdge=0;iEdge<4;iEdge++){
+ int p1 = indices[iEdge];
+ int p2 = indices[(iEdge+1)%4];
+ int nbP = order-1;
+ if (p1 == 0 && p2 == 1){for (int i=4+0*nbP;i< 4+1*nbP;i++) tmp[index++] = vtcs[i];}
+ else if (p1 == 1 && p2 == 2){for (int i=4+1*nbP;i< 4+2*nbP;i++) tmp[index++] = vtcs[i];}
+ else if (p1 == 2 && p2 == 3){for (int i=4+2*nbP;i< 4+3*nbP;i++) tmp[index++] = vtcs[i];}
+ else if (p1 == 3 && p2 == 0){for (int i=4+3*nbP;i< 4+4*nbP;i++) tmp[index++] = vtcs[i];}
+
+ else if (p1 == 1 && p2 == 0){for (int i=4+1*nbP-1;i>= 4+0*nbP;i--) tmp[index++] = vtcs[i];}
+ else if (p1 == 2 && p2 == 1){for (int i=4+2*nbP-1;i>= 4+1*nbP;i--) tmp[index++] = vtcs[i];}
+ else if (p1 == 3 && p2 == 2){for (int i=4+3*nbP-1;i>= 4+2*nbP;i--) tmp[index++] = vtcs[i];}
+ else if (p1 == 0 && p2 == 3){for (int i=4+4*nbP-1;i>= 4+3*nbP;i--) tmp[index++] = vtcs[i];}
+ else printf("ouuls\n");
+ }
+ for (int i=0;i<index;i++)vtcs[start+4+i] = tmp[i];
+ start += (4+index);
+ }
+
+ order -= 2;
+ if (start >= vtcs.size()) break;
+ }
+ // printf("after : ");
+ // for (int i=0;i<vtcs.size();i++)printf("%p ",vtcs[i]);
+ // printf("\n");
+}
+
+
// KH: check face orientation wrt element ...
static void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &vf,
@@ -506,53 +571,100 @@ static void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &v
{
fullMatrix<double> points;
int start = 0;
- switch (nPts){
- case 0 :
- case 1 :
- // do nothing (e.g. for 2nd order tri faces or for quad faces)
- break;
- case 2 :
- points = polynomialBases::find(MSH_TRI_10)->points;
- start = 9;
- break;
- case 3 :
- points = polynomialBases::find(MSH_TRI_15)->points;
- start = 12;
- break;
- case 4 :
- points = polynomialBases::find(MSH_TRI_21)->points;
- start = 15;
- break;
- case 5 :
- points = polynomialBases::find(MSH_TRI_28)->points;
- start = 18;
- break;
- case 6 :
- points = polynomialBases::find(MSH_TRI_36)->points;
- start = 21;
- break;
- case 7 :
- points = polynomialBases::find(MSH_TRI_45)->points;
- start = 24;
- break;
- case 8 :
- points = polynomialBases::find(MSH_TRI_55)->points;
- start = 27;
- break;
- case 9 :
- points = polynomialBases::find(MSH_TRI_66)->points;
- start = 30;
- break;
- default :
- Msg::Error("getFaceVertices not implemented for order %i\n",nPts +1);
- break;
+
+ switch(ele->getType()){
+ case TYPE_TET:
+ {
+ switch (nPts){
+ case 0 :
+ case 1 :
+ // do nothing (e.g. for 2nd order tri faces)
+ break;
+ case 2 :
+ points = polynomialBases::find(MSH_TRI_10)->points;
+ start = 9;
+ break;
+ case 3 :
+ points = polynomialBases::find(MSH_TRI_15)->points;
+ start = 12;
+ break;
+ case 4 :
+ points = polynomialBases::find(MSH_TRI_21)->points;
+ start = 15;
+ break;
+ case 5 :
+ points = polynomialBases::find(MSH_TRI_28)->points;
+ start = 18;
+ break;
+ case 6 :
+ points = polynomialBases::find(MSH_TRI_36)->points;
+ start = 21;
+ break;
+ case 7 :
+ points = polynomialBases::find(MSH_TRI_45)->points;
+ start = 24;
+ break;
+ case 8 :
+ points = polynomialBases::find(MSH_TRI_55)->points;
+ start = 27;
+ break;
+ case 9 :
+ points = polynomialBases::find(MSH_TRI_66)->points;
+ start = 30;
+ break;
+ default :
+ Msg::Error("getFaceVertices not implemented for order %i\n",nPts +1);
+ break;
+ }
+ break;
+ }
+ case TYPE_HEX:
+ {
+ switch (nPts){
+ case 0 :
+ // do nothing (e.g. for 1nd order quad faces)
+ break;
+ case 1 :
+ points = polynomialBases::find(MSH_QUA_9)->points;
+ break;
+ case 2 :
+ points = polynomialBases::find(MSH_QUA_16)->points;
+ break;
+ case 3 :
+ points = polynomialBases::find(MSH_QUA_25)->points;
+ break;
+ case 4 :
+ points = polynomialBases::find(MSH_QUA_36)->points;
+ break;
+ case 5 :
+ points = polynomialBases::find(MSH_QUA_49)->points;
+ break;
+ case 6 :
+ points = polynomialBases::find(MSH_QUA_64)->points;
+ break;
+ case 7 :
+ points = polynomialBases::find(MSH_QUA_81)->points;
+ break;
+ case 8 :
+ points = polynomialBases::find(MSH_QUA_100)->points;
+ break;
+ default :
+ Msg::Error("getFaceVertices not implemented for order %i\n",nPts +1);
+ break;
+ }
+ start = (nPts+2)*(nPts+2) - nPts*nPts;
+ // printf("start = %d\n",start);
+ break;
+ }
}
-
+
for(int i = 0; i < ele->getNumFaces(); i++){
MFace face = ele->getFace(i);
+ // printf(" face %d %d vertices\n",i,face.getNumVertices());
faceContainer::iterator fIter = faceVertices.find(face);
- if (fIter != faceVertices.end()) {
+ if (fIter != faceVertices.end() ) {
std::vector<MVertex*> vtcs = fIter->second;
+ // printf("found %d\n",vtcs.size());
if(face.getNumVertices() == 3 && nPts > 1){ // tri face
int orientation;
bool swap;
@@ -562,7 +674,16 @@ static void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &v
Msg::Error("Error in face lookup for recuperation of high order face nodes");
}
else if(face.getNumVertices() == 4){ // quad face
- // TODO reorient if more than 1 face vertex
+ int orientation;
+ bool swap;
+ if (fIter->first.computeCorrespondence(face, orientation, swap)){
+ // printf("%d %d %d %d\n",face.getVertex(0)->getNum(),face.getVertex(1)->getNum(),face.getVertex(2)->getNum(),face.getVertex(3)->getNum());
+ // printf("%d %d %d %d\n",fIter->first.getVertex(0)->getNum(),fIter->first.getVertex(1)->getNum(),fIter->first.getVertex(2)->getNum(),fIter->first.getVertex(3)->getNum());
+ // printf("%d %d\n",orientation,swap);
+ reorientQuadPoints(vtcs, orientation, swap, nPts-1);
+ }
+ else
+ Msg::Error("Error in face lookup for recuperation of high order face nodes");
}
vf.insert(vf.end(), vtcs.begin(), vtcs.end());
}
@@ -600,17 +721,17 @@ static void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &v
}
}
else if(face.getNumVertices() == 4){ // quad face
- for(int j = 0; j < nPts; j++){
- for(int k = 0; k < nPts; k++){
+ // printf(" quad face %d pts\n",nPts);
+ // FIXME : CURVED Quad Face
+ for (int k = start; k < points.size1(); k++) {
// parameters are between -1 and 1
- double t1 = 2. * (double)(j + 1) / (nPts + 1) - 1.;
- double t2 = 2. * (double)(k + 1) / (nPts + 1) - 1.;
- SPoint3 pc = face.interpolate(t1, t2);
- MVertex *v = new MVertex(pc.x(), pc.y(), pc.z(), gr);
- vtcs.push_back(v);
- gr->mesh_vertices.push_back(v);
- vf.push_back(v);
- }
+ double t1 = points(k, 0);
+ double t2 = points(k, 1);
+ SPoint3 pc = face.interpolate(t1, t2);
+ MVertex *v = new MVertex(pc.x(), pc.y(), pc.z(), gr);
+ vtcs.push_back(v);
+ gr->mesh_vertices.push_back(v);
+ vf.push_back(v);
}
}
}
@@ -623,37 +744,81 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
fullMatrix<double> points;
int start = 0;
- switch (nPts){
- case 0:
- case 1:
- case 2:
- // done: return!
- return;
- case 3 :
- points = polynomialBases::find(MSH_TET_35)->points;
- break;
- case 4 :
- points = polynomialBases::find(MSH_TET_56)->points;
- break;
- case 5 :
- points = polynomialBases::find(MSH_TET_84)->points;
- break;
- case 6 :
- points = polynomialBases::find(MSH_TET_120)->points;
- break;
- case 7 :
- points = polynomialBases::find(MSH_TET_165)->points;
- break;
- case 8 :
- points = polynomialBases::find(MSH_TET_220)->points;
- break;
- case 9 :
- points = polynomialBases::find(MSH_TET_286)->points;
- break;
- default :
- Msg::Error("getRegionVertices not implemented for order %i\n", nPts+1);
+ switch (incomplete->getType()){
+ case TYPE_TET :
+ {
+ switch (nPts){
+ case 0:
+ case 1:
+ return;
+ case 2:
+ points = polynomialBases::find(MSH_TET_20)->points;
+ break;
+ case 3 :
+ points = polynomialBases::find(MSH_TET_35)->points;
+ break;
+ case 4 :
+ points = polynomialBases::find(MSH_TET_56)->points;
+ break;
+ case 5 :
+ points = polynomialBases::find(MSH_TET_84)->points;
+ break;
+ case 6 :
+ points = polynomialBases::find(MSH_TET_120)->points;
+ break;
+ case 7 :
+ points = polynomialBases::find(MSH_TET_165)->points;
+ break;
+ case 8 :
+ points = polynomialBases::find(MSH_TET_220)->points;
+ break;
+ case 9 :
+ points = polynomialBases::find(MSH_TET_286)->points;
+ break;
+ default :
+ Msg::Error("getRegionVertices not implemented for order %i\n", nPts+1);
+ }
+ start = ((nPts+2)*(nPts+3)*(nPts+4)-(nPts-2)*(nPts-1)*(nPts))/6;
+ break;
+ }
+ case TYPE_HEX :
+ {
+ switch (nPts){
+ case 0:
+ return;
+ case 1:
+ points = polynomialBases::find(MSH_HEX_27)->points;
+ break;
+ case 2 :
+ points = polynomialBases::find(MSH_HEX_64)->points;
+ break;
+ case 3 :
+ points = polynomialBases::find(MSH_HEX_125)->points;
+ break;
+ case 4 :
+ points = polynomialBases::find(MSH_HEX_216)->points;
+ break;
+ case 5 :
+ points = polynomialBases::find(MSH_HEX_343)->points;
+ break;
+ case 6 :
+ points = polynomialBases::find(MSH_HEX_512)->points;
+ break;
+ case 7 :
+ points = polynomialBases::find(MSH_HEX_729)->points;
+ break;
+ case 8 :
+ points = polynomialBases::find(MSH_HEX_1000)->points;
+ break;
+ default :
+ Msg::Error("getRegionVertices not implemented for order %i\n", nPts+1);
+ }
+ start = (nPts+2)*(nPts+2)*(nPts+2) - (nPts)*(nPts)*(nPts) ;
+ break;
+ }
}
- start = ((nPts+2)*(nPts+3)*(nPts+4)-(nPts-2)*(nPts-1)*(nPts))/6;
+
+ // printf ("incomplete type %d starts at %d\n",incomplete->getTypeForMSH(),start);
for(int k = start; k < points.size1(); k++){
MVertex *v;
@@ -662,6 +827,7 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
const double t3 = points(k, 2);
SPoint3 pos;
incomplete->pnt(t1,t2,t3,pos);
+ // printf("pnt(%g %g %g) = %g %g %g\n",t1,t2,t3,pos.x(),pos.y(),pos.z());
v = new MVertex(pos.x(),pos.y(),pos.z(),gr);
gr->mesh_vertices.push_back(v);
vr.push_back(v);
@@ -841,27 +1007,69 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
std::vector<MHexahedron*> hexahedra2;
for(unsigned int i = 0; i < gr->hexahedra.size(); i++){
MHexahedron *h = gr->hexahedra[i];
- std::vector<MVertex*> ve, vf;
+ h->revert();
+ std::vector<MVertex*> ve, vf, vr;
getEdgeVertices(gr, h, ve, edgeVertices, linear, nPts, displ2D, displ3D);
- if(incomplete){
- hexahedra2.push_back
- (new MHexahedron20(h->getVertex(0), h->getVertex(1), h->getVertex(2),
- h->getVertex(3), h->getVertex(4), h->getVertex(5),
- h->getVertex(6), h->getVertex(7), ve[0], ve[1], ve[2],
- ve[3], ve[4], ve[5], ve[6], ve[7], ve[8], ve[9], ve[10],
- ve[11]));
+ if(nPts == 1){
+ if(incomplete){
+ hexahedra2.push_back
+ (new MHexahedron20(h->getVertex(0), h->getVertex(1), h->getVertex(2),
+ h->getVertex(3), h->getVertex(4), h->getVertex(5),
+ h->getVertex(6), h->getVertex(7), ve[0], ve[1], ve[2],
+ ve[3], ve[4], ve[5], ve[6], ve[7], ve[8], ve[9], ve[10],
+ ve[11]));
+ }
+ else{
+ getFaceVertices(gr, h, vf, faceVertices, edgeVertices, linear, nPts);
+ SPoint3 pc = h->barycenter();
+ MVertex *v = new MVertex(pc.x(), pc.y(), pc.z(), gr);
+ gr->mesh_vertices.push_back(v);
+ hexahedra2.push_back
+ (new MHexahedron27(h->getVertex(0), h->getVertex(1), h->getVertex(2),
+ h->getVertex(3), h->getVertex(4), h->getVertex(5),
+ h->getVertex(6), h->getVertex(7), ve[0], ve[1], ve[2],
+ ve[3], ve[4], ve[5], ve[6], ve[7], ve[8], ve[9], ve[10],
+ ve[11], vf[0], vf[1], vf[2], vf[3], vf[4], vf[5], v));
+ }
}
- else{
+ else {
+ // printf("coucou %d\n",ve.size());
getFaceVertices(gr, h, vf, faceVertices, edgeVertices, linear, nPts);
- SPoint3 pc = h->barycenter();
- MVertex *v = new MVertex(pc.x(), pc.y(), pc.z(), gr);
- gr->mesh_vertices.push_back(v);
- hexahedra2.push_back
- (new MHexahedron27(h->getVertex(0), h->getVertex(1), h->getVertex(2),
- h->getVertex(3), h->getVertex(4), h->getVertex(5),
- h->getVertex(6), h->getVertex(7), ve[0], ve[1], ve[2],
- ve[3], ve[4], ve[5], ve[6], ve[7], ve[8], ve[9], ve[10],
- ve[11], vf[0], vf[1], vf[2], vf[3], vf[4], vf[5], v));
+ ve.insert(ve.end(), vf.begin(), vf.end());
+ // printf("coucou %d\n",ve.size());
+ MHexahedronN incpl(h->getVertex(0), h->getVertex(1), h->getVertex(2), h->getVertex(3),
+ h->getVertex(4), h->getVertex(5), h->getVertex(6), h->getVertex(7),
+ ve, nPts + 1);
+ getRegionVertices(gr, &incpl, h, vr, linear, nPts);
+
+ // printf("oulax\n");
+
+ ve.insert(ve.end(), vr.begin(), vr.end());
+ // printf("coucou %d\n",ve.size());
+ MHexahedron* n = new MHexahedronN(h->getVertex(0), h->getVertex(1), h->getVertex(2), h->getVertex(3),
+ h->getVertex(4), h->getVertex(5), h->getVertex(6), h->getVertex(7),
+ ve, nPts + 1);
+
+ /*
+ FILE *F = fopen ("etst.pos","w");
+ fprintf(F,"View \"\"{\n");
+ for (int i=0;i<n->getNumVertices();i++){
+ // printf("%3d %12.5E %12.5E %12.5E\n",n->getVertex(i)->getNum(),n->getVertex(i)->x(),n->getVertex(i)->y(),n->getVertex(i)->z());
+ fprintf(F,"SP(%12.5E,%12.5E,%12.5E){%3d};\n",n->getVertex(i)->x(),n->getVertex(i)->y(),n->getVertex(i)->z(),i);//n->getVertex(i)->getNum());
+ }
+ fprintf(F,"};\n");
+ fclose(F);
+ */
+ /*
+ const polynomialBasis* temp = n->getFunctionSpace();
+ for (int i=0;i<n->getNumVertices();i++){
+ SPoint3 pt;
+ incpl.pnt(temp->points(i,0),temp->points(i,1),temp->points(i,2),pt);
+ printf("%12.5E %12.5E %12.5E -- %12.5E %12.5E %12.5E -- %12.5E %12.5E %12.5E\n",n->getVertex(i)->x(),n->getVertex(i)->y(),n->getVertex(i)->z(),
+ pt.x(),pt.y(),pt.z(),pt.x()-n->getVertex(i)->x(),pt.y()-n->getVertex(i)->y(),pt.z()-n->getVertex(i)->z());
+ }
+ */
+ hexahedra2.push_back(n);
}
delete h;
}
@@ -1113,9 +1321,11 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
int nPts = order - 1;
- Msg::StatusBar(2, true, "Meshing order %d...", order);
+ if (!linear)Msg::StatusBar(2, true, "Meshing order %d, curvilinear ON...", order);
+ else Msg::StatusBar(2, true, "Meshing order %d, curvilinear OFF...", order);
double t1 = Cpu();
+
// first, make sure to remove any existsing second order vertices/elements
SetOrder1(m);
@@ -1143,6 +1353,7 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
// we do that model face by model face
std::vector<MElement*> bad;
double worst;
+ printJacobians(m, "smoothness_b.pos");
if (displ2D){
checkHighOrderTriangles("Before optimization", m, bad, worst);
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
@@ -1166,7 +1377,7 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
if(displ2D) delete displ2D;
if(displ3D) delete displ3D;
- // printJacobians(m, "smoothness.pos");
+ printJacobians(m, "smoothness.pos");
double t2 = Cpu();
Msg::StatusBar(2, true, "Done meshing order %d (%g s)", order, t2 - t1);
diff --git a/Mesh/highOrderSmoother.cpp b/Mesh/highOrderSmoother.cpp
index 9dea04cb76c51824508ca0c4bf0b761356aff14c..5afed29d69d9ac1ae414947e56adafbda16dcf70 100644
--- a/Mesh/highOrderSmoother.cpp
+++ b/Mesh/highOrderSmoother.cpp
@@ -29,6 +29,7 @@
#include "elasticityTerm.h"
#include "linearSystemCSR.h"
#include "linearSystemFull.h"
+#include "linearSystemPetsc.h"
#define SQU(a) ((a)*(a))
@@ -323,6 +324,28 @@ void addOneLayer(const std::vector<MElement*> &v,
}
}
+
+void _printJacobiansAtNodes (const char * name, std::vector<MElement*> &v)
+{
+ FILE *fd = fopen (name,"w");
+ fprintf(fd,"$MeshFormat\n2 0 8\n$EndMeshFormat\n$ElementNodeData\n1\n\"det J\"\n1\n0.0\n3\n1\n1\n%d\n",v.size());
+ for (int i=0; i<v.size(); i++){
+ const polynomialBasis* pb = v[i]->getFunctionSpace();
+ fprintf(fd,"%d %d",v[i]->getNum(),v[i]->getNumVertices());
+ for (int j=0; j < v[i]->getNumVertices(); j++){
+ const double _u = pb->points(j,0);
+ const double _v = pb->points(j,1);
+ const double _w = pb->points(j,2);
+ double jac[3][3];
+ double J = v[i]->getJacobian (_u,_v,_w, jac);
+ fprintf(fd," %g",J);
+ }
+ fprintf(fd,"\n");
+ }
+ fprintf(fd,"$EndElementNodeData\n");
+ fclose(fd);
+}
+
void highOrderSmoother::smooth(GFace *gf, bool metric)
{
std::vector<MElement*> v;
@@ -331,6 +354,14 @@ void highOrderSmoother::smooth(GFace *gf, bool metric)
v.insert(v.begin(), gf->quadrangles.begin(),gf->quadrangles.end());
Msg::Info("Smoothing high order mesh : model face %d (%d elements)", gf->tag(),
v.size());
+
+ gf->model()->writeMSH("before_smoothing.msh");
+ _printJacobiansAtNodes ("before_smoothing_jac.msh", v);
+ smooth_with_mixed_formulation(v,0.2);
+ gf->model()->writeMSH("after_smoothing.msh");
+ _printJacobiansAtNodes ("after_smoothing_jac.msh", v);
+ return;
+
if (metric)smooth_metric(v,gf);
else smooth(v);
}
@@ -620,6 +651,196 @@ double highOrderSmoother::smooth_metric_(std::vector<MElement*> & v,
return dx;
}
+
+double highOrderSmoother::apply_incremental_displacement (double max_incr,
+ std::vector<MElement*> & v,
+ bool mixed,
+ double thres,
+ char *meshName,
+ std::vector<MElement*> & disto)
+{
+#ifdef HAVE_PETSC
+ // assume that the mesh is OK, yet already curved
+ linearSystemPETSc<double> *lsys = new linearSystemPETSc<double>;
+ dofManager<double> myAssembler(lsys);
+ elasticityMixedTerm El_mixed (0, 1.0, .333, getTag());
+ elasticityTerm El (0, 1.0, .333, getTag());
+
+ std::set<MVertex*> _vertices;
+
+ //+++++++++ Boundary Conditions & Numbering +++++++++++++++++++++++++++++++
+ // fix all dof that correspond to vertices on the boundary
+ // the value is equal
+ for (unsigned int i = 0; i < v.size(); i++){
+ for (int j = 0; j < v[i]->getNumVertices(); j++){
+ MVertex *vert = v[i]->getVertex(j);
+ _vertices.insert(vert);
+ }
+ }
+
+ //+++++++++ Fix d tr(eps) = 0 +++++++++++++++++++++++++++++++
+ if (mixed){
+ for (unsigned int i = 0; i < disto.size(); i++){
+ for (int j = 0; j < disto[i]->getNumVertices(); j++){
+ MVertex *vert = disto[i]->getVertex(j);
+ myAssembler.fixVertex(vert, 4, getTag(), 0.0);
+ }
+ }
+ }
+
+ for (std::set<MVertex*>::iterator it = _vertices.begin(); it != _vertices.end(); ++it){
+ MVertex *vert = *it;
+ std::map<MVertex*,SVector3>::iterator itt = _targetLocation.find(vert);
+ // impose displacement @ boundary
+ if (vert->onWhat()->geomType() != GEntity::Line && vert->onWhat()->tag() != 1){ // TEST
+ if (itt != _targetLocation.end() && vert->onWhat()->dim() < _dim){
+ myAssembler.fixVertex(vert, 0, getTag(), itt->second.x()-vert->x());
+ myAssembler.fixVertex(vert, 1, getTag(), itt->second.y()-vert->y());
+ myAssembler.fixVertex(vert, 2, getTag(), itt->second.z()-vert->z());
+ }
+ // ensure we do not touch any vertex that is on the boundary
+ else if (vert->onWhat()->dim() < _dim){
+ myAssembler.fixVertex(vert, 0, getTag(), 0);
+ myAssembler.fixVertex(vert, 1, getTag(), 0);
+ myAssembler.fixVertex(vert, 2, getTag(), 0);
+ }
+ }
+ if (_dim == 2)myAssembler.fixVertex(vert, 2, getTag(), 0);
+ // number vertices
+ myAssembler.numberVertex(vert, 0, getTag());
+ myAssembler.numberVertex(vert, 1, getTag());
+ myAssembler.numberVertex(vert, 2, getTag());
+ if (mixed){
+ myAssembler.numberVertex(vert, 3, getTag());
+ myAssembler.numberVertex(vert, 4, getTag());
+ }
+ }
+
+ //+++++++++ Assembly & Solve ++++++++++++++++++++++++++++++++++++
+ if (myAssembler.sizeOfR()){
+ // assembly of the elasticity term on the
+ for (unsigned int i = 0; i < v.size(); i++){
+ SElement se(v[i]);
+ if (mixed)El_mixed.addToMatrix(myAssembler, &se);
+ else El.addToMatrix(myAssembler, &se);
+ }
+ // solve the system
+ lsys->systemSolve();
+ }
+
+ //+++++++++ Move vertices @ maximum ++++++++++++++++++++++++++++++++++++
+ FILE *fd = fopen ("d.msh","w");
+ fprintf(fd,"$MeshFormat\n2 0 8\n$EndMeshFormat\n$NodeData\n1\n\"tr(sigma)\"\n1\n0.0\n3\n1\n3\n%d\n",_vertices.size());
+ for (std::set<MVertex*>::iterator it = _vertices.begin(); it != _vertices.end(); ++it){
+ double ax, ay, az;
+ myAssembler.getDofValue(*it, 0, getTag(), ax);
+ myAssembler.getDofValue(*it, 1, getTag(), ay);
+ myAssembler.getDofValue(*it, 2, getTag(), az);
+ (*it)->x() += max_incr*ax;
+ (*it)->y() += max_incr*ay;
+ (*it)->z() += max_incr*az;
+ fprintf(fd,"%d %g %g %g\n",(*it)->getIndex(), ax,ay,az);
+ }
+ fprintf(fd,"$EndNodeData\n");
+ fclose(fd);
+
+ //+------------------- Check now if elements are ok ---------------+++++++
+
+ (*_vertices.begin())->onWhat()->model()->writeMSH(meshName);
+
+ /*
+ std::map<MVertex*,double> facts;
+ while(1){
+ int count = 0;
+ for (unsigned int i = 0; i < v.size(); i++){
+ double disto = v[i]->distoShapeMeasure();
+ if (disto < thres){
+ count ++;
+ for (int j = 0; j < v[i]->getNumVertices(); j++){
+ MVertex *vert = v[i]->getVertex(j);
+ std::map<MVertex*,double>::iterator it = facts.find(vert);
+ if (it == facts.end())facts[vert] = .1;
+ else if (it->second < 1){
+ it->second += .1;
+ double ax, ay, az;
+ myAssembler.getDofValue(vert, 0, getTag(), ax);
+ myAssembler.getDofValue(vert, 1, getTag(), ay);
+ myAssembler.getDofValue(vert, 2, getTag(), az);
+ vert->x() -= .1*ax;
+ vert->y() -= .1*ay;
+ vert->z() -= .1*az;
+ }
+ }
+ }
+ if (!count)break;
+ Msg::Info("%d elements are negative : reducing displacement",count);
+ }
+ */
+
+ double percentage = max_incr * 100.;
+ while(1){
+ std::vector<MElement*> disto;
+ double minD;
+ getDistordedElements(v, 0.5, disto, minD);
+ if (minD < thres){
+ percentage -= 10.;
+ for (std::set<MVertex*>::iterator it = _vertices.begin(); it != _vertices.end(); ++it){
+ double ax, ay, az;
+ myAssembler.getDofValue(*it, 0, getTag(), ax);
+ myAssembler.getDofValue(*it, 1, getTag(), ay);
+ myAssembler.getDofValue(*it, 2, getTag(), az);
+ (*it)->x() -= .1*ax;
+ (*it)->y() -= .1*ay;
+ (*it)->z() -= .1*az;
+ }
+ }
+ else break;
+ }
+
+ delete lsys;
+ return percentage;
+#endif
+ return 0.0;
+}
+
+int highOrderSmoother::smooth_with_mixed_formulation (std::vector<MElement*> &all,
+ double alpha)
+{
+ int ITER = 0;
+ double minD, FACT = 1.0;
+ std::vector<MElement*> disto;
+ // move the points to their straight sided locations
+ for (unsigned int i = 0; i < all.size(); i++)
+ moveToStraightSidedLocation(all[i]);
+ // apply the displacement
+ double percentage = 0.0;
+ while(1){
+ char NN[256];
+ sprintf(NN,"smoothing-%d.msh",ITER++);
+ //double percentage_of_what_is_left = apply_incremental_displacement (1.,all, false, alpha,NN,disto);
+ double percentage_of_what_is_left = apply_incremental_displacement (1.,all, true, alpha,NN,all);
+ percentage += (1.-percentage) * percentage_of_what_is_left/100.;
+ Msg::Info("The smoother was able to do %3d percent of the motion",(int)(percentage*100.));
+ if (percentage_of_what_is_left == 0.0) break;
+ else if (percentage_of_what_is_left == 100.)break;
+ }
+
+ getDistordedElements(all, 0.3, disto, minD);
+ Msg::Info("iter %d : %d elements / %d distorted min Disto = %g ",ITER,
+ all.size(), disto.size(), minD);
+
+ if (0 && percentage < 0.99){
+ char NN[256];
+ sprintf(NN,"smoothing-%d.msh",ITER++);
+ double percentage_of_what_is_left = apply_incremental_displacement (1.0,all, true, alpha,NN,disto);
+ percentage += (1.-percentage) * percentage_of_what_is_left/100.;
+ Msg::Info("The mixed smoother was able to do %3d percent of the motion",(int)(percentage*100.));
+ }
+
+ return 1;
+}
+
+
void highOrderSmoother::smooth(std::vector<MElement*> &all)
{
#ifdef HAVE_TAUCS
@@ -641,7 +862,7 @@ void highOrderSmoother::smooth(std::vector<MElement*> &all)
Msg::Info("%d elements / %d distorted min Disto = %g\n",
all.size(), v.size(), minD);
- if (!v.size()) return;
+ if (!v.size());
const int nbLayers = 6;
for (int i = 0; i < nbLayers; i++){
diff --git a/Mesh/highOrderSmoother.h b/Mesh/highOrderSmoother.h
index 9ff58bc77043b8d448697d07bfcffe0454472219..26a684cc08fc396d35d346f85841964fe51b0eef 100644
--- a/Mesh/highOrderSmoother.h
+++ b/Mesh/highOrderSmoother.h
@@ -44,6 +44,11 @@ public:
_straightSidedLocation[v] = d;
_targetLocation[v] = SPoint3(v->x(),v->y(),v->z());
}
+ int smooth_with_mixed_formulation(std::vector<MElement*> & ,
+ double alpha);
+ double apply_incremental_displacement (double max_incr, std::vector<MElement*> & v,
+ bool mixed, double thres, char *meshName,
+ std::vector<MElement*> & disto);
void smooth(std::vector<MElement*> & );
double smooth_metric_(std::vector<MElement*> &, GFace *gf,
dofManager<double> &myAssembler,
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index 07956971eb07e8882fafb64d2654effd3f78c934..1b5c8942fe142df3d8b42987ba499841cfa72a83 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -626,7 +626,7 @@ static int _quadWithOneVertexOnBoundary (GFace *gf,
MVertex *v2,
MVertex *v3,
MVertex *v4){
- return 0;
+ //return 0;
if (v1->onWhat()->dim() < 2 &&
v2->onWhat()->dim() == 2 &&
v3->onWhat()->dim() == 2 &&
@@ -1862,10 +1862,10 @@ void recombineIntoQuads(GFace *gf,
if(gf->geomType() == GEntity::DiscreteSurface && !gf->getCompound())
haveParam = false;
- gf->model()->writeMSH("before.msh");
if(haveParam && topologicalOpti)
removeFourTrianglesNodes(gf, false);
+ gf->model()->writeMSH("before.msh");
int success = _recombineIntoQuads(gf, 0);
// gf->addLayersOfQuads(1, 0);
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index 5a9738360601c81093024b37a871cd9214f7c349..2b2e4dc88f984cfd2e972dc151998154b4c365aa 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -235,6 +235,7 @@ double qmDistorsionOfMapping (MTriangle *e)
const double u = pts[i].pt[0];
const double v = pts[i].pt[1];
const double di = mesh_functional_distorsion (e, u, v);
+ // printf("di = %g\n",di);
dmin = (i == 0)? di : std::min(dmin, di);
}
const fullMatrix<double>& points = e->getFunctionSpace()->points;
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 81be8ed044ecfe4b3f402b49d1187a0a78d5c302..59b59416dfe2634417f6dfae31a3aa7b26767e4e 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -159,7 +159,7 @@ const fullMatrix<double> &polynomialBasis::getGradientAtFaceIntegrationPoints(in
const polynomialBasis *fsFace = polynomialBases::find(closures[iClosure].type);
gaussIntegration::get(fsFace->parentType, integrationOrder, integrationFace, weight);
fullMatrix<double> integration(integrationFace.size1(), 3);
- double f[256];
+ double f[1256];
for (int i = 0; i < integrationFace.size1(); i++){
fsFace->f(integrationFace(i,0),integrationFace(i,1),integrationFace(i,2),f);
for(size_t j=0; j<closures[iClosure].size(); j++) {
@@ -170,7 +170,7 @@ const fullMatrix<double> &polynomialBasis::getGradientAtFaceIntegrationPoints(in
}
fullMatrix<double> &v = dfAtFace[iClosure];
v.resize(nbPsi, integration.size1()*3);
- double g[256][3];
+ double g[1256][3];
for (size_t xi = 0; xi< integration.size1(); xi++) {
df(integration(xi,0 ), integration(xi,1), integration(xi,2), g);
for (int j = 0; j < nbPsi; j++) {
@@ -203,6 +203,7 @@ static fullMatrix<double> generatePascalQuad(int order)
return monomials;
}
+
/*
00 10 20 30 40 ⋯
01 11 21 31 41 ⋯
@@ -213,26 +214,25 @@ static fullMatrix<double> generatePascalQuad(int order)
*/
// generate all monomials xi^m * eta^n with n and m <= order
-static fullMatrix<double> generatePascalHex(int order)
+static fullMatrix<double> generatePascalHex(int order, bool serendip)
{
-
- fullMatrix<double> monomials( (order+1)*(order+1)*(order+1), 3);
+ int siz = (order+1)*(order+1)*(order+1);
+ if (serendip) siz -= (order-1)*(order-1)*(order-1);
+ fullMatrix<double> monomials( siz, 3);
int index = 0;
- for (int p = 0; p <= order; p++) {
- for(int i = 0; i < p; i++) {
- for(int j = 0; j < i; j++, index++) {
- monomials(index, 0) = p;
- monomials(index, 1) = i;
- monomials(index, 2) = j;
- }
- }
- for(int i = 0; i <= p; i++, index++) {
- for(int j = 0; j <= p; i++, index++) {
- monomials(index, 0) = p-i;
- monomials(index, 1) = p;
+ for (int i = 0; i <= order; i++) {
+ for (int j = 0; j <= order; j++) {
+ for (int k = 0; k <= order; k++) {
+ if (!serendip || i<2 || j<2 || k<2){
+ monomials(index,0) = i;
+ monomials(index,1) = j;
+ monomials(index,2) = k;
+ index++;
+ }
}
}
}
+ // printf("siz %d %d\n",siz,index );
return monomials;
}
@@ -265,6 +265,7 @@ static fullMatrix<double> generatePascalQuadSerendip(int order)
return monomials;
}
+
/*static fullMatrix<double> generatePascalQuadSerendip(int order)
{
@@ -843,7 +844,8 @@ static fullMatrix<double> gmshGeneratePointsQuad(int order, bool serendip)
point(index, 1) = point(p0, 1) + i*(point(p1,1)-point(p0,1))/order;
}
}
- if (order > 2 && !serendip) {
+ // FIXIT it was > 2 and I beleive it is >= 2 !!
+ if (order >= 2 && !serendip) {
fullMatrix<double> inner = gmshGeneratePointsQuad(order - 2, false);
inner.scale(1. - 2./order);
point.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
@@ -857,6 +859,120 @@ static fullMatrix<double> gmshGeneratePointsQuad(int order, bool serendip)
return point;
}
+static fullMatrix<double> gmshGeneratePointsHex(int order, bool serendip)
+{
+ int nbPoints = (order+1)*(order+1)*(order+1);
+ if (serendip) nbPoints -= (order-1)*(order-1)*(order-1);
+ fullMatrix<double> point(nbPoints, 3);
+
+ // should be a public member of MHexahedron, just copied
+ static const int edges[12][2] = {
+ {0, 1},
+ {0, 3},
+ {0, 4},
+ {1, 2},
+ {1, 5},
+ {2, 3},
+ {2, 6},
+ {3, 7},
+ {4, 5},
+ {4, 7},
+ {5, 6},
+ {6, 7}
+ };
+ static const int faces[6][4] = {
+ {0, 3, 2, 1},
+ {0, 1, 5, 4},
+ {0, 4, 7, 3},
+ {1, 2, 6, 5},
+ {2, 3, 7, 6},
+ {4, 5, 6, 7}
+ };
+
+ if (order > 0) {
+ point(0, 0) = -1;
+ point(0, 1) = -1;
+ point(0, 2) = -1;
+ point(1, 0) = 1;
+ point(1, 1) = -1;
+ point(1, 2) = -1;
+ point(2, 0) = 1;
+ point(2, 1) = 1;
+ point(2, 2) = -1;
+ point(3, 0) = -1;
+ point(3, 1) = 1;
+ point(3, 2) = -1;
+
+ point(4, 0) = -1;
+ point(4, 1) = -1;
+ point(4, 2) = 1;
+ point(5, 0) = 1;
+ point(5, 1) = -1;
+ point(5, 2) = 1;
+ point(6, 0) = 1;
+ point(6, 1) = 1;
+ point(6, 2) = 1;
+ point(7, 0) = -1;
+ point(7, 1) = 1;
+ point(7, 2) = 1;
+
+ if (order > 1) {
+ int index = 8;
+ for (int iedge=0; iedge<12; iedge++) {
+ int p0 = edges[iedge][0];
+ int p1 = edges[iedge][1];
+ for (int i = 1; i < order; i++, index++) {
+ point(index, 0) = point(p0, 0) + i*(point(p1,0)-point(p0,0))/order;
+ point(index, 1) = point(p0, 1) + i*(point(p1,1)-point(p0,1))/order;
+ point(index, 2) = point(p0, 2) + i*(point(p1,2)-point(p0,2))/order;
+ }
+ }
+
+ fullMatrix<double> fp = gmshGeneratePointsQuad(order - 2, false);
+ fp.scale(1. - 2./order);
+ for (int iface=0; iface<6; iface++) {
+ int p0 = faces[iface][0];
+ int p1 = faces[iface][1];
+ int p2 = faces[iface][2];
+ int p3 = faces[iface][3];
+ for (int i=0;i<fp.size1();i++, index++){
+ const double u = fp(i,0);
+ const double v = fp(i,1);
+ const double newU =
+ 0.25*(1.-u)*(1.-v)*point(p0,0) +
+ 0.25*(1.+u)*(1.-v)*point(p1,0) +
+ 0.25*(1.+u)*(1.+v)*point(p2,0) +
+ 0.25*(1.-u)*(1.+v)*point(p3,0) ;
+ const double newV =
+ 0.25*(1.-u)*(1.-v)*point(p0,1) +
+ 0.25*(1.+u)*(1.-v)*point(p1,1) +
+ 0.25*(1.+u)*(1.+v)*point(p2,1) +
+ 0.25*(1.-u)*(1.+v)*point(p3,1) ;
+ const double newW =
+ 0.25*(1.-u)*(1.-v)*point(p0,2) +
+ 0.25*(1.+u)*(1.-v)*point(p1,2) +
+ 0.25*(1.+u)*(1.+v)*point(p2,2) +
+ 0.25*(1.-u)*(1.+v)*point(p3,2) ;
+ point(index, 0) = newU;
+ point(index, 1) = newV;
+ point(index, 2) = newW;
+ }
+ }
+ if (!serendip) {
+ fullMatrix<double> inner = gmshGeneratePointsHex(order - 2, false);
+ inner.scale(1. - 2./order);
+ point.copy(inner, 0, nbPoints - index, 0, 3, index, 0);
+ }
+ }
+ }
+ else if (order == 0){
+ point(0, 0) = 0;
+ point(0, 1) = 0;
+ point(0, 2) = 0;
+ }
+ return point;
+}
+
static fullMatrix<double> generateLagrangeMonomialCoefficients
(const fullMatrix<double>& monomial, const fullMatrix<double>& point)
{
@@ -1273,8 +1389,17 @@ static void generateClosureOrder0(polynomialBasis::clCont &closure, int nb)
}
}
-std::map<int, polynomialBasis> polynomialBases::fs;
+void _oulala_ (const char* a, fullMatrix<double> & m){
+ FILE *f = fopen (a,"w");
+ fprintf(f,"View \"\"{\n");
+ for (int i=0;i<m.size1();i++){
+ fprintf(f,"SP(%g,%g,%g){%d};\n",m(i,0),m(i,1),m(i,2),i);
+ }
+ fprintf(f,"};\n");
+ fclose (f);
+}
+std::map<int, polynomialBasis> polynomialBases::fs;
const polynomialBasis *polynomialBases::find(int tag)
{
std::map<int, polynomialBasis>::const_iterator it = fs.find(tag);
@@ -1860,6 +1985,142 @@ const polynomialBasis *polynomialBases::find(int tag)
generateFaceClosurePrism(F.closures, 2);
generateFaceClosurePrismFull(F.fullClosures, F.closureRef, 2);
break;
+ case MSH_HEX_20 :
+ F.numFaces = 6;
+ F.monomials = generatePascalHex(2,true);
+ F.points = gmshGeneratePointsHex(2, true);
+ F.parentType = TYPE_HEX;
+ // generateFaceClosureHex(F.closures, 2);
+ // generateFaceClosureHexFull(F.fullClosures, F.closureRef, 2);
+ break;
+ case MSH_HEX_27 :
+ F.numFaces = 6;
+ F.monomials = generatePascalHex(2,false);
+ F.points = gmshGeneratePointsHex(2, false);
+ F.parentType = TYPE_HEX;
+ // generateFaceClosureHex(F.closures, 2);
+ // generateFaceClosureHexFull(F.fullClosures, F.closureRef, 2);
+ break;
+ case MSH_HEX_56 :
+ // printf("in complete hex at order 3\n");
+ F.numFaces = 6;
+ F.monomials = generatePascalHex(3,true);
+ F.points = gmshGeneratePointsHex(3, true);
+ F.parentType = TYPE_HEX;
+ // _oulala_("hex56.pos",F.points);
+ // generateFaceClosureHex(F.closures, 2);
+ // generateFaceClosureHexFull(F.fullClosures, F.closureRef, 2);
+ break;
+ case MSH_HEX_64 :
+ // printf("complete hex at order 3\n");
+ F.numFaces = 6;
+ F.monomials = generatePascalHex(3,false);
+ // printf("pas complete hex at order 3 done\n");
+ F.points = gmshGeneratePointsHex(3, false);
+ // printf("complete hex at order 3 done\n");
+ // F.points.print("P3");
+ F.parentType = TYPE_HEX;
+ // _oulala_("hex64.pos",F.points);
+ // generateFaceClosureHex(F.closures, 2);
+ // generateFaceClosureHexFull(F.fullClosures, F.closureRef, 2);
+ break;
+ case MSH_HEX_98 :
+ F.numFaces = 6;
+ F.monomials = generatePascalHex(4,true);
+ F.points = gmshGeneratePointsHex(4, true);
+ F.parentType = TYPE_HEX;
+ // generateFaceClosureHex(F.closures, 2);
+ // generateFaceClosureHexFull(F.fullClosures, F.closureRef, 2);
+ break;
+ case MSH_HEX_125 :
+ F.numFaces = 6;
+ F.monomials = generatePascalHex(4,false);
+ F.points = gmshGeneratePointsHex(4, false);
+ F.parentType = TYPE_HEX;
+ // _oulala_("hex125.pos",F.points);
+ // generateFaceClosureHex(F.closures, 2);
+ // generateFaceClosureHexFull(F.fullClosures, F.closureRef, 2);
+ break;
+ case MSH_HEX_152 :
+ F.numFaces = 6;
+ F.monomials = generatePascalHex(5,true);
+ F.points = gmshGeneratePointsHex(5, true);
+ F.parentType = TYPE_HEX;
+ // generateFaceClosureHex(F.closures, 2);
+ // generateFaceClosureHexFull(F.fullClosures, F.closureRef, 2);
+ break;
+ case MSH_HEX_216 :
+ F.numFaces = 6;
+ F.monomials = generatePascalHex(5,false);
+ F.points = gmshGeneratePointsHex(5, false);
+ F.parentType = TYPE_HEX;
+ // generateFaceClosureHex(F.closures, 2);
+ // generateFaceClosureHexFull(F.fullClosures, F.closureRef, 2);
+ break;
+ case MSH_HEX_222 :
+ F.numFaces = 6;
+ F.monomials = generatePascalHex(6,true);
+ F.points = gmshGeneratePointsHex(6, true);
+ F.parentType = TYPE_HEX;
+ // generateFaceClosureHex(F.closures, 2);
+ // generateFaceClosureHexFull(F.fullClosures, F.closureRef, 2);
+ break;
+ case MSH_HEX_343 :
+ F.numFaces = 6;
+ F.monomials = generatePascalHex(6,false);
+ F.points = gmshGeneratePointsHex(6, false);
+ F.parentType = TYPE_HEX;
+ // generateFaceClosureHex(F.closures, 2);
+ // generateFaceClosureHexFull(F.fullClosures, F.closureRef, 2);
+ break;
+ case MSH_HEX_296 :
+ F.numFaces = 6;
+ F.monomials = generatePascalHex(7,true);
+ F.points = gmshGeneratePointsHex(7, true);
+ F.parentType = TYPE_HEX;
+ // generateFaceClosureHex(F.closures, 2);
+ // generateFaceClosureHexFull(F.fullClosures, F.closureRef, 2);
+ break;
+ case MSH_HEX_512 :
+ F.numFaces = 6;
+ F.monomials = generatePascalHex(7,false);
+ F.points = gmshGeneratePointsHex(7, false);
+ F.parentType = TYPE_HEX;
+ // generateFaceClosureHex(F.closures, 2);
+ // generateFaceClosureHexFull(F.fullClosures, F.closureRef, 2);
+ break;
+ case MSH_HEX_386 :
+ F.numFaces = 6;
+ F.monomials = generatePascalHex(8,true);
+ F.points = gmshGeneratePointsHex(8, true);
+ F.parentType = TYPE_HEX;
+ // generateFaceClosureHex(F.closures, 2);
+ // generateFaceClosureHexFull(F.fullClosures, F.closureRef, 2);
+ break;
+ case MSH_HEX_729 :
+ F.numFaces = 6;
+ F.monomials = generatePascalHex(8,false);
+ F.points = gmshGeneratePointsHex(8, false);
+ F.parentType = TYPE_HEX;
+ // generateFaceClosureHex(F.closures, 2);
+ // generateFaceClosureHexFull(F.fullClosures, F.closureRef, 2);
+ break;
+ case MSH_HEX_488 :
+ F.numFaces = 6;
+ F.monomials = generatePascalHex(9,true);
+ F.points = gmshGeneratePointsHex(9, true);
+ F.parentType = TYPE_HEX;
+ // generateFaceClosureHex(F.closures, 2);
+ // generateFaceClosureHexFull(F.fullClosures, F.closureRef, 2);
+ break;
+ case MSH_HEX_1000 :
+ F.numFaces = 6;
+ F.monomials = generatePascalHex(9,false);
+ F.points = gmshGeneratePointsHex(9, false);
+ F.parentType = TYPE_HEX;
+ // generateFaceClosureHex(F.closures, 2);
+ // generateFaceClosureHexFull(F.fullClosures, F.closureRef, 2);
+ break;
default :
Msg::Error("Unknown function space %d: reverting to TET_4", tag);
F.numFaces = 4;
@@ -1900,7 +2161,7 @@ const fullMatrix<double> &polynomialBases::findInjector(int tag1, int tag2)
fullMatrix<double> inj(fs1.points.size1(), fs2.points.size1());
- double sf[256];
+ double sf[1256];
for (int i = 0; i < fs1.points.size1(); i++) {
fs2.f(fs1.points(i, 0), fs1.points(i, 1), fs1.points(i, 2), sf);
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index c1481612e0693ac42ecc342cf671016d8fd186ed..0170b2c63306bbf605aa6905d73bf9d5eb21af56 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -115,7 +115,7 @@ class polynomialBasis
}
inline void f(double u, double v, double w, double *sf) const
{
- double p[256];
+ double p[1256];
evaluateMonomials(u, v, w, p);
for (int i = 0; i < coefficients.size1(); i++) {
sf[i] = 0;
@@ -128,7 +128,7 @@ class polynomialBasis
// implemented) and is easier to bind
inline void f(fullMatrix<double> &coord, fullMatrix<double> &sf)
{
- double p[256];
+ double p[1256];
sf.resize (coord.size1(), coefficients.size1());
for (int iPoint=0; iPoint< coord.size1(); iPoint++) {
evaluateMonomials(coord(iPoint,0), coord(iPoint,1), coord(iPoint,2), p);
diff --git a/Solver/elasticityTerm.cpp b/Solver/elasticityTerm.cpp
index 7c67875577213c20956dcb7728455f77b5a3da36..f03ec416ac293709a08b6b263f74b75889e27b73 100644
--- a/Solver/elasticityTerm.cpp
+++ b/Solver/elasticityTerm.cpp
@@ -29,10 +29,10 @@ void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
double C12 = FACT * _nu / (1 - 2 * _nu);
double C44 = (C11 - C12) / 2;
// FIXME : PLANE STRESS !!!
- FACT = _E / (1-_nu*_nu);
- C11 = FACT;
- C12 = _nu * FACT;
- C44 = (1.-_nu)*.5*FACT;
+ //FACT = _E / (1-_nu*_nu);
+ //C11 = FACT;
+ //C12 = _nu * FACT;
+ //C44 = (1.-_nu)*.5*FACT;
const double C[6][6] =
{ {C11, C12, C12, 0, 0, 0},
@@ -130,3 +130,168 @@ void elasticityTerm::elementVector(SElement *se, fullVector<double> &m) const
}
}
}
+
+/*
+
+B_d is the deviatoric part of the FE strains
+
+K_{uu} = \int_\Omega B^T_d H B_d dv
+
+*/
+
+void elasticityMixedTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
+{
+ setPolynomialBasis(se);
+ MElement *e = se->getMeshElement();
+ int integrationOrder = 3 * (e->getPolynomialOrder()) + 2;
+ int npts;
+ IntPt *GP;
+ double jac[3][3];
+ double invjac[3][3];
+ double Grads[100][3];
+ double SF[100];
+ e->getIntegrationPoints(integrationOrder, &npts, &GP);
+
+ m.setAll(0.);
+ fullMatrix<double> KUU(3 * _sizeN,3 * _sizeN);
+ fullMatrix<double> KUP(3 * _sizeN,_sizeM);
+ fullMatrix<double> KUG(3 * _sizeN,_sizeM);
+ fullMatrix<double> KPG(_sizeM,_sizeM);
+ fullMatrix<double> KGG(_sizeM,_sizeM);
+
+ double FACT = _E / ((1 + _nu)*(1.-2*_nu));
+ double C11 = FACT * (1 - _nu) ;
+ double C12 = FACT * _nu ;
+ double C44 = FACT* (1 - 2.*_nu)*.5;
+ double _K = 3*_E / (1 - 2 * _nu);
+ // FIXME : PLANE STRESS !!!
+ // FACT = _E / (1-_nu*_nu);
+ // C11 = FACT;
+ // C12 = _nu * FACT;
+ // C44 = (1.-_nu)*.5*FACT;
+
+ const double C[6][6] =
+ { {C11, C12, C12, 0, 0, 0},
+ {C12, C11, C12, 0, 0, 0},
+ {C12, C12, C11, 0, 0, 0},
+ { 0, 0, 0, C44, 0, 0},
+ { 0, 0, 0, 0, C44, 0},
+ { 0, 0, 0, 0, 0, C44} };
+
+ fullMatrix<double> H (6, 6);
+ fullMatrix<double> B (6, 3 * _sizeN);
+ fullMatrix<double> BDEV (6, 3 * _sizeN);
+ fullMatrix<double> BDIL (1, 3 * _sizeN);
+ fullMatrix<double> BDILT (3 * _sizeN,1);
+ fullMatrix<double> BTH(3 * _sizeN, 6);
+ fullMatrix<double> BT(3 * _sizeN, 6);
+ fullMatrix<double> trBTH(3 * _sizeN, 1);
+ fullMatrix<double> MT(_sizeM, 1);
+ fullMatrix<double> M(1,_sizeM);
+ for (int i = 0; i < 6; i++)
+ for (int j = 0; j < 6; j++)
+ H(i, j) = C[i][j];
+
+ double _dimension = 3.0;
+
+ for (int i = 0; i < npts; i++){
+ const double u = GP[i].pt[0];
+ const double v = GP[i].pt[1];
+ const double w = GP[i].pt[2];
+ const double weight = GP[i].weight;
+ const double detJ = e->getJacobian(u, v, w, jac);
+ inv3x3(jac, invjac);
+ se->gradNodalShapeFunctions(u, v, w, invjac,Grads);
+ e->getShapeFunctions (u,v,w,SF,_polyOrderM);
+
+ B.setAll(0.);
+ BT.setAll(0.);
+ for (int j = 0; j < _sizeN; j++){
+
+ // printf("%g %g %g\n",Grads[j][0],Grads[j][1],Grads[j][2]);
+
+ BDILT(j, 0) = BDIL(0, j) = Grads[j][0]/_dimension;
+
+ BT(j, 0) = B(0, j) = Grads[j][0]-Grads[j][0]/_dimension;
+ BT(j, 1) = B(1, j) = -Grads[j][1]/_dimension;
+ BT(j, 2) = B(2, j) = -Grads[j][2]/_dimension;
+
+ BT(j, 3) = B(3, j) = Grads[j][1];
+ BT(j, 4) = B(4, j) = Grads[j][2];
+
+ BDILT(j + _sizeN, 0) = BDIL(0, j + _sizeN) = Grads[j][1]/_dimension;
+
+ BT(j + _sizeN, 0) = B(0, j + _sizeN) = -Grads[j][0]/_dimension;
+ BT(j + _sizeN, 1) = B(1, j + _sizeN) = Grads[j][1]-Grads[j][1]/_dimension;
+ BT(j + _sizeN, 2) = B(2, j + _sizeN) = -Grads[j][2]/_dimension;
+
+ BT(j + _sizeN, 3) = B(3, j + _sizeN) = Grads[j][0];
+ BT(j + _sizeN, 5) = B(5, j + _sizeN) = Grads[j][2];
+
+ BDILT(j + 2 * _sizeN, 0) = BDIL(0, j + 2 * _sizeN) = Grads[j][2]/_dimension;
+
+ BT(j + 2 * _sizeN, 0) = B(0, j + 2 * _sizeN) = -Grads[j][0]/_dimension;
+ BT(j + 2 * _sizeN, 1) = B(1, j + 2 * _sizeN) = -Grads[j][1]/_dimension;
+ BT(j + 2 * _sizeN, 2) = B(2, j + 2 * _sizeN) = Grads[j][2]-Grads[j][2]/_dimension;
+
+ BT(j + 2 * _sizeN, 4) = B(4, j + 2 * _sizeN) = Grads[j][1];
+ BT(j + 2 * _sizeN, 5) = B(5, j + 2 * _sizeN) = Grads[j][0];
+ }
+
+ for (int j = 0; j < _sizeM; j++){
+ MT(j,0) = M(0,j) = SF[j];
+ }
+
+ // printf("BT (%d %d) x H (%d,%d) = BTH(%d,%d)\n",BT.size1(),BT.size2(), H.size1(),H.size2(),BTH.size1(),BTH.size2());
+
+ BTH.setAll(0.);
+ BTH.gemm(BT, H);
+
+ for (int j = 0; j < 3*_sizeN; j++){
+ trBTH(j,0) = BTH(j,0) + BTH(j,1) + BTH(j,2) ;
+ }
+
+ KUU.gemm(BTH, B, weight * detJ);
+ KUP.gemm(BDILT, M, weight * detJ);
+ KPG.gemm(MT, M, -weight * detJ);
+ KUG.gemm(trBTH, M, weight * detJ/_dimension);
+ KGG.gemm(MT, M, weight * detJ * (_K)/(_dimension*_dimension));
+
+ }
+ /*
+ 3*_sizeN _sizeM _sizeM
+
+ KUU KUP KUG
+ m = KUP^t 0 KPG
+ KUG^T KPG^t KGG
+
+ */
+
+ for (int i=0;i<3*_sizeN;i++){
+ for (int j=0;j<3*_sizeN;j++)
+ m(i,j) = KUU(i,j); // assemble KUU
+
+ for (int j=0;j<_sizeM;j++){
+ m(i,j+3*_sizeN) = KUP(i,j); // assemble KUP
+ m(j+3*_sizeN,i) = KUP(i,j); // assemble KUP^T
+ }
+
+ for (int j=0;j<_sizeM;j++){
+ m(i,j+3*_sizeN+_sizeM) = KUG(i,j); // assemble KUG
+ m(j+3*_sizeN+_sizeM,i) = KUG(i,j); // assemble KUG^t
+ }
+ }
+ for (int i=0;i<_sizeM;i++){
+ for (int j=0;j<_sizeM;j++){
+ m(i+3*_sizeN,j+3*_sizeN+_sizeM) = KPG(i,j); // assemble KPG
+ m(j+3*_sizeN+_sizeM,i+3*_sizeN) = KPG(i,j); // assemble KPG^t
+
+ m(i+3*_sizeN+_sizeM,j+3*_sizeN+_sizeM) = KGG(i,j); // assemble KGG
+
+ }
+ }
+ // m.print("Mixed");
+ return;
+
+}
+
diff --git a/Solver/elasticityTerm.h b/Solver/elasticityTerm.h
index c89b1593d749d85c91d2f9f10dfc84af84786968..e167fd4a61a65f0261b939741548f1651a1e4952 100644
--- a/Solver/elasticityTerm.h
+++ b/Solver/elasticityTerm.h
@@ -9,6 +9,7 @@
#include "femTerm.h"
#include "Gmsh.h"
#include "GModel.h"
+#include "polynomialBasis.h"
#include "SElement.h"
#include "fullMatrix.h"
@@ -57,4 +58,74 @@ class elasticityTerm : public femTerm<double> {
void elementVector(SElement *se, fullVector<double> &m) const;
};
+/*
+ Formulation of elasticity with 3 fields
+ -) displacement U
+ -) lagrange multipliers p and g
+
+ g = trace (epsilon)
+ p = trace (epsilon)
+
+
+*/
+
+class elasticityMixedTerm : public femTerm<double> {
+protected:
+ double _E, _nu;
+ mutable int _iField, _polyOrderN, _polyOrderM, _sizeN, _sizeM;
+ mutable polynomialBasis *_pN,*_pM;
+ void setPolynomialBasis (SElement *se) const{
+ _polyOrderN = se->getMeshElement()->getPolynomialOrder();
+ _polyOrderM = se->getMeshElement()->getPolynomialOrder();
+ _pN = (polynomialBasis*)se->getMeshElement()->getFunctionSpace(_polyOrderN);
+ _pM = (polynomialBasis*)se->getMeshElement()->getFunctionSpace(_polyOrderM);
+ _sizeN = _pN->coefficients.size1();
+ _sizeM = _pM->coefficients.size1();
+ }
+public:
+ void setField(int i){_iField = i;}
+ // element matrix size : 3 dofs per vertex
+ virtual int sizeOfR(SElement *se) const
+ {
+ setPolynomialBasis(se);
+ return 3 * _sizeN + 2 * _sizeM;
+ }
+ virtual int sizeOfC(SElement *se) const
+ {
+ return sizeOfR(se);
+ }
+ // order dofs in the local matrix :
+ // dx1, dx2 ... dxn, dy1, ..., dyn, dz1, ... dzn ,
+ // p1,p2,...,pm, g1,g2,...,gm
+ Dof getLocalDofR(SElement *se, int iRow) const
+ {
+ setPolynomialBasis(se);
+ MElement *e = se->getMeshElement();
+ int iComp;
+ int ithLocalVertex;
+ if (iRow < 3 * _sizeN){
+ iComp = iRow / _sizeN;
+ ithLocalVertex = iRow % _sizeN;
+ }
+ else {
+ iRow -= 3 * _sizeN;
+ iComp = 3 + iRow / _sizeM;
+ ithLocalVertex = iRow % _sizeM;
+ }
+ return Dof(e->getShapeFunctionNode(ithLocalVertex)->getNum(),
+ Dof::createTypeWithTwoInts(iComp, _iField));
+ }
+ Dof getLocalDofC(SElement *se, int iCol) const
+ {
+ return getLocalDofR(se,iCol);
+ }
+ public:
+ elasticityMixedTerm(GModel *gm, double E, double nu, int field)
+ : femTerm<double>(gm), _E(E), _nu(nu), _iField(field) {}
+ void elementMatrix(SElement *se, fullMatrix<double> &m) const;
+ void elementVector(SElement *se, fullVector<double> &m) const {}
+ void setYoung (double E){_E = E;}
+};
+
+
#endif
diff --git a/Solver/functionSpace.h b/Solver/functionSpace.h
index 644378f0270a6e555612e4c1efcd0c3dd2a170f6..ec0da235b6a2e8d73a8dde19cb4859ed06dc51d0 100644
--- a/Solver/functionSpace.h
+++ b/Solver/functionSpace.h
@@ -237,7 +237,7 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
if(ele->getParent()) ele = ele->getParent();
int ndofs= ele->getNumShapeFunctions();
vals.reserve(vals.size()+ndofs);
- double valsuvw[256];
+ double valsuvw[1256];
ele->getShapeFunctions(u, v, w, valsuvw);
for(int i = 0; i < ndofs; ++i)
vals.push_back(valsuvw[i]);
diff --git a/benchmarks/3d/Sphere.geo b/benchmarks/3d/Sphere.geo
index b48c3ef31aa6b4df8b2d2a607a3592326e1aaf51..3e445e48d0de1595b4ceaf6b65f12a87afb1e4ed 100644
--- a/benchmarks/3d/Sphere.geo
+++ b/benchmarks/3d/Sphere.geo
@@ -1,6 +1,6 @@
// Gmsh project created on Fri Feb 1 12:17:44 2008
-Delete All;
+//Delete All;
RSphere = 1.;
lcSphere = .25;
diff --git a/benchmarks/extrude/u_shape_boundary_layer.geo b/benchmarks/extrude/u_shape_boundary_layer.geo
index 1c157cf1f155a6733d349570445bcd412032dc24..6b66c16d38f84d78ef171d42519badc7332c2ba1 100644
--- a/benchmarks/extrude/u_shape_boundary_layer.geo
+++ b/benchmarks/extrude/u_shape_boundary_layer.geo
@@ -11,8 +11,8 @@ BSpline(1) = {5, 4, 3, 2, 1};
BSpline(2) = {1, 9, 8, 7, 5};
Line(3) = {5, 6};
-Extrude { Line{1,-3}; Layers{5,0.2}; Using Index[0]; }
-Extrude { Line{2,3}; Layers{5,0.2}; Using Index[1]; }
+Extrude { Line{1,-3}; Layers{5,0.04}; Using Index[0]; }
+Extrude { Line{2,3}; Layers{5,0.04}; Using Index[1]; }
// fix leading edge by hand
Coherence Point {25, 16};
@@ -27,4 +27,4 @@ Line(22) = {34, 31};
Line(23) = {31, 32};
Line Loop(24) = {22, 23, 20, 21};
Line Loop(25) = {4, 12, 16, -18, 9, 8};
-Plane Surface(26) = {24, 25};
+//Plane Surface(26) = {24, 25};