diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index d729b27519de6534c979485e02c653215f27e904..14243e6d6c84bd432fbb18f75cf51fcecefe37d8 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1839,17 +1839,6 @@ GModel *GModel::buildCutGModel(gLevelset *ls, bool cutElem, bool saveTri)
Msg::Info("Mesh cutting completed (%g s)", Cpu() - t1);
- //emi debug
- // std::vector<GEntity*> entities;
- // cutGM->getEntities(entities);
- // for(unsigned int i = 0; i < entities.size(); i++){
- // std::vector<int> phys = entities[i]->getPhysicalEntities();
- // for (int j= 0; j < phys.size(); j++){
- // std::string name = cutGM->getPhysicalName(entities[i]->dim(), phys[j]);
- // printf("dim =%d elem=%d phys =%s \n", entities[i]->dim(),entities[i]->tag(), name.c_str() );
- // }
- // }
-
return cutGM;
}
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index fe21cb7cd93c30948848157bd1d624020ae71618..d9ac6d2db2f41a4e0a0b082878b5521a4ce606e6 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -49,7 +49,7 @@ const polynomialBasis* MQuadrangle::getFunctionSpace(int o) const
case 8: return polynomialBases::find(MSH_QUA_81);
case 9: return polynomialBases::find(MSH_QUA_100);
case 10: return polynomialBases::find(MSH_QUA_121);
- default: Msg::Error("Order %d triangle function space not implemented", order);
+ default: Msg::Error("Order %d quadrangle function space not implemented", order);
}
return 0;
}
@@ -85,7 +85,7 @@ const JacobianBasis* MQuadrangle::getJacobianFuncSpace(int o) const
case 8: return JacobianBasis::find(MSH_QUA_81);
case 9: return JacobianBasis::find(MSH_QUA_100);
case 10: return JacobianBasis::find(MSH_QUA_121);
- default: Msg::Error("Order %d triangle function space not implemented", order);
+ default: Msg::Error("Order %d quadrangle function space not implemented", order);
}
return 0;
}
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 18b2caf491d553d135483d7afdc16629baa8c3fc..28374f7c0b81c1ec6c18394033b8cd8a617ee38e 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -1399,7 +1399,7 @@ const polynomialBasis *polynomialBases::find(int tag)
case MSH_LIN_1 : F.parentType = TYPE_LIN; F.order = 0; F.serendip = false; break;
case MSH_LIN_2 : F.parentType = TYPE_LIN; F.order = 1; F.serendip = false; break;
case MSH_LIN_3 : F.parentType = TYPE_LIN; F.order = 2; F.serendip = false; break;
- case MSH_LIN_4 : F.parentType = TYPE_LIN; F.order = 3; F.serendip = false; break;
+ case MSH_LIN_4 : printf("line 4 \n"); F.parentType = TYPE_LIN; F.order = 3; F.serendip = false; break;
case MSH_LIN_5 : F.parentType = TYPE_LIN; F.order = 4; F.serendip = false; break;
case MSH_LIN_6 : F.parentType = TYPE_LIN; F.order = 5; F.serendip = false; break;
case MSH_LIN_7 : F.parentType = TYPE_LIN; F.order = 6; F.serendip = false; break;
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index 0110c5f087efdc4e68740a0e2e531dd7c481d772..87a64b4c1a5ea6c716e22ae5321a3927ca20e744 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -93,6 +93,7 @@ class polynomialBasis
fullMatrix<double> monomials;
fullMatrix<double> coefficients;
int numFaces;
+ inline int getNumShapeFunctions() const {return coefficients.size1();};
// for a given face/edge, with both a sign and a rotation, give an
// ordered list of nodes on this face/edge
inline const std::vector<int> &getClosure(int id) const
@@ -127,7 +128,7 @@ class polynomialBasis
double p[1256];
evaluateMonomials(u, v, w, p);
for (int i = 0; i < coefficients.size1(); i++) {
- sf[i] = 0;
+ sf[i] = 0.0;
for (int j = 0; j < coefficients.size2(); j++) {
sf[i] += coefficients(i, j) * p[j];
}
diff --git a/Solver/function.cpp b/Solver/function.cpp
index d65b84c3c207a29e36ac1f98cd29dd278ef25de8..b30c9aa9b68d1b9198b6471d6a624e23f576610d 100644
--- a/Solver/function.cpp
+++ b/Solver/function.cpp
@@ -497,7 +497,7 @@ class functionLevelsetSmooth : public function {
}
functionLevelsetSmooth(const function *f0, const double valMin, const double valPlus, const double E) : function(f0->getNbCol())
{
- printf("Levelset bandwidth is E = %g \n", E);
+ //printf("Levelset bandwidth is E = %g \n", E);
setArgument (_f0, f0);
_valMin = valMin;
_valPlus = valPlus;
diff --git a/contrib/DiscreteIntegration/Integration3D.cpp b/contrib/DiscreteIntegration/Integration3D.cpp
index 76fb9977c67e1675b7cdb8523d9e631ab77ac995..4a31c0254bf57ab9ebcb13ab99cfc23666977aca 100644
--- a/contrib/DiscreteIntegration/Integration3D.cpp
+++ b/contrib/DiscreteIntegration/Integration3D.cpp
@@ -4,6 +4,8 @@
#include "Integration3D.h"
#include "recurCut.h"
#include "Gauss.h"
+#include "polynomialBasis.h"
+#include "GmshDefines.h"
#define ZERO_LS_TOL 1.e-3
#define EQUALITY_TOL 1.e-15
@@ -613,9 +615,24 @@ DI_Element & DI_Element::operator= (const DI_Element &rhs){
}
return *this;
}
+void DI_Element::getShapeFunctions(double u, double v, double w, double s[], int o) const
+{
+ //printf("type elem =%d order =%d\n", type(), o);
+ const polynomialBasis* fs = getFunctionSpace(o);
+ if(fs) fs->f(u, v, w, s);
+ else Msg::Error("Function space not implemented for this type of element");
+}
+
+void DI_Element::getGradShapeFunctions(double u, double v, double w, double s[][3],
+ int o) const
+{
+ const polynomialBasis* fs = getFunctionSpace(o);
+ if(fs) fs->df(u, v, w, s);
+ else Msg::Error("Function space not implemented for this type of element");
+}
void DI_Element::setPolynomialOrder (int o) {
if(polOrder_ == o) return;
- delete [] mid_;
+ if (mid_) delete [] mid_;
polOrder_ = o;
switch (o) {
case 1 :
@@ -632,12 +649,13 @@ void DI_Element::setPolynomialOrder (int o) {
}
return;
default :
- printf("Order %d line function space not implemented ", o); print();
+ return;
+ //printf("Order %d line function space not implemented ", o); print();
}
}
void DI_Element::setPolynomialOrder (int o, const DI_Element *e, const std::vector<gLevelset *> &RPNi) {
if(polOrder_ == o) return;
- delete [] mid_;
+ if (mid_) delete [] mid_;
polOrder_ = o;
switch (o) {
case 1 :
@@ -742,6 +760,7 @@ bool DI_Element::addQuadEdge (int edge, DI_Point *xm,
printf("detJ<0 when trying to add a quadratic edge in "); print();
return false;
}
+ printf("in add quad edge \n");
computeIntegral();
return true;
}
@@ -1042,370 +1061,160 @@ bool DI_ElementLessThan::operator()(const DI_Element *e1, const DI_Element *e2)
}
// DI_Line methods --------------------------------------------------------------------------------
-void DI_Line::computeIntegral() {
- switch (getPolynomialOrder()) {
- case 1 :
- integral_ = LineLength(pt(0), pt(1));
- break;
- case 2 :
- integral_ = LineLength(pt(0), mid(0)) + LineLength(mid(0), pt(1));
- break;
- default :
- printf("Order %d line function space not implemented ", getPolynomialOrder()); print();
- }
-}
-void DI_Line::getShapeFunctions (double u, double v, double w, double s[], int ord) const {
- int order = (ord == -1) ? getPolynomialOrder() : ord;
- double valm = (fabs(1. - u) < 1.e-16) ? 0. : (1. - u);
- double valp = (fabs(1. + u) < 1.e-16) ? 0. : (1. + u);
+const polynomialBasis* DI_Line::getFunctionSpace(int o) const{
+ int order = (o == -1) ? getPolynomialOrder() : o;
switch (order) {
- case 1 :
- s[0] = valm / 2.;
- s[1] = valp / 2.;
- break;
- case 2 :
- s[0] = -u * valm / 2.;
- s[1] = u * valp / 2.;
- s[2] = valm * valp;
- break;
- default : printf("Order %d line function space not implemented ", order); print();
- }
-}
-void DI_Line::getGradShapeFunctions (const double u, const double v, const double w, double s[][3], int ord) const {
- int order = (ord == -1) ? getPolynomialOrder() : ord;
- switch (order) {
- case 1 :
- s[0][0] = -0.5; s[0][1] = 0.; s[0][2] = 0.;
- s[1][0] = 0.5; s[1][1] = 0.; s[1][2] = 0.;
- break;
- case 2 :
- s[0][0] = u - 0.5; s[0][1] = 0.; s[0][2] = 0.;
- s[1][0] = u + 0.5; s[1][1] = 0.; s[1][2] = 0.;
- s[2][0] = -2. * u; s[2][1] = 0.; s[2][2] = 0.;
- break;
- default :
- printf("Order %d line function space not implemented ", order); print();
+ case 0: return polynomialBases::find(MSH_LIN_1);
+ case 1: return polynomialBases::find(MSH_LIN_2);
+ case 2: return polynomialBases::find(MSH_LIN_3);
+ case 3: return polynomialBases::find(MSH_LIN_4);
+ case 4: return polynomialBases::find(MSH_LIN_5);
+ case 5: return polynomialBases::find(MSH_LIN_6);
+ case 6: return polynomialBases::find(MSH_LIN_7);
+ case 7: return polynomialBases::find(MSH_LIN_8);
+ case 8: return polynomialBases::find(MSH_LIN_9);
+ case 9: return polynomialBases::find(MSH_LIN_10);
+ case 10: return polynomialBases::find(MSH_LIN_11);
+ default: Msg::Error("Order %d line function space not implemented", order);
}
+ return 0;
}
-// DI_Triangle methods ----------------------------------------------------------------------------
-void DI_Triangle::computeIntegral() {
- switch (getPolynomialOrder()) {
- case 1 :
- integral_ = TriSurf(pt(0), pt(1), pt(2));
- break;
- case 2 :
- integral_ = TriSurf(pt(0), mid(0), mid(2)) + TriSurf(pt(1), mid(0), mid(1))
- + TriSurf(pt(2), mid(1), mid(2)) + TriSurf(mid(0), mid(1), mid(2));
- break;
- default :
- printf("Order %d triangle function space not implemented ", getPolynomialOrder()); print();
- }
+void DI_Line::computeIntegral() {
+ integral_ = LineLength(pt(0), pt(1));
+ // switch (getPolynomialOrder()) {
+ // case 1 :
+ // integral_ = LineLength(pt(0), pt(1));
+ // break;
+ // case 2 :
+ // integral_ = LineLength(pt(0), mid(0)) + LineLength(mid(0), pt(1));
+ // break;
+ // default :
+ // printf("Order %d line function space not implemented ", getPolynomialOrder()); print();
+ // }
}
-void DI_Triangle::getShapeFunctions (double u, double v, double w, double s[], int ord) const {
- int order = (ord == -1) ? getPolynomialOrder() : ord;
- double val1 = (fabs(1. - u - v) < 1.e-16) ? 0. : (1. - u - v);
+// DI_Triangle methods ----------------------------------------------------------------------------
+const polynomialBasis* DI_Triangle::getFunctionSpace(int o) const
+{
+ int order = (o == -1) ? getPolynomialOrder() : o;
switch (order) {
- case 1 :
- s[0] = val1;
- s[1] = u;
- s[2] = v;
- break;
- case 2 :
- s[0] = val1 * (1. - 2. * u - 2. * v);
- s[1] = u * (2. * u - 1.);
- s[2] = v * (2. * v - 1.);
- s[3] = 4. * u * val1;
- s[4] = 4. * u * v;
- s[5] = 4. * v * val1;
- break;
- default :
- printf("Order %d triangle function space not implemented ", order); print();
- }
+ case 0: return polynomialBases::find(MSH_TRI_1);
+ case 1: return polynomialBases::find(MSH_TRI_3);
+ case 2: return polynomialBases::find(MSH_TRI_6);
+ case 3: return polynomialBases::find(MSH_TRI_10);
+ case 4: return polynomialBases::find(MSH_TRI_15);
+ case 5: return polynomialBases::find(MSH_TRI_21);
+ case 6: return polynomialBases::find(MSH_TRI_28);
+ case 7: return polynomialBases::find(MSH_TRI_36);
+ case 8: return polynomialBases::find(MSH_TRI_45);
+ case 9: return polynomialBases::find(MSH_TRI_55);
+ case 10: return polynomialBases::find(MSH_TRI_66);
+ default: Msg::Error("Order %d triangle function space not implemented", order);
+ }
}
-void DI_Triangle::getGradShapeFunctions (const double u, const double v, const double w,
- double s[][3], int ord) const {
- int order = (ord == -1) ? getPolynomialOrder() : ord;
- switch (order) {
- case 1 :
- s[0][0] = -1.; s[0][1] = -1.; s[0][2] = 0.;
- s[1][0] = 1.; s[1][1] = 0.; s[1][2] = 0.;
- s[2][0] = 0.; s[2][1] = 1.; s[2][2] = 0.;
- break;
- case 2 :
- s[0][0] = 4. * u + 4. * v - 3.;
- s[0][1] = 4. * u + 4. * v - 3.;
- s[0][2] = 0.;
- s[1][0] = 4. * u -1.;
- s[1][1] = 0.;
- s[1][2] = 0.;
- s[2][0] = 0.;
- s[2][1] = 4. * v - 1.;
- s[2][2] = 0.;
- s[3][0] = -8. * u - 4. * v + 4.;
- s[3][1] = -4. * u;
- s[3][2] = 0.;
- s[4][0] = 4. * v;
- s[4][1] = 4. * u;
- s[4][2] = 0.;
- s[5][0] = -4. * v;
- s[5][1] = -4. * u - 8. * v + 4.;
- s[5][2] = 0.;
- break;
- default :
- printf("Order %d triangle function space not implemented ", order); print();
- }
+void DI_Triangle::computeIntegral() {
+ integral_ = TriSurf(pt(0), pt(1), pt(2));
+ // switch (getPolynomialOrder()) {
+ // case 1 :
+ // integral_ = TriSurf(pt(0), pt(1), pt(2));
+ // break;
+ // case 2 :
+ // integral_ = TriSurf(pt(0), mid(0), mid(2)) + TriSurf(pt(1), mid(0), mid(1))
+ // + TriSurf(pt(2), mid(1), mid(2)) + TriSurf(mid(0), mid(1), mid(2));
+ // break;
+ // default :
+ // printf("Order %d triangle function space not implemented ", getPolynomialOrder()); print();
+ // }
}
double DI_Triangle::quality() const {
return qualityTri(pt(0), pt(1), pt(2));
}
// DI_Quad methods --------------------------------------------------------------------------------
-void DI_Quad::computeIntegral() {
- switch (getPolynomialOrder()) {
- case 1 :
- integral_ = TriSurf(pt(0), pt(1), pt(2)) + TriSurf(pt(0), pt(2), pt(3));
- break;
- case 2 :
- integral_ = TriSurf(pt(0), mid(0), mid(4)) + TriSurf(pt(0), mid(4), mid(3))
- + TriSurf(pt(1), mid(1), mid(4)) + TriSurf(pt(1), mid(4), mid(0))
- + TriSurf(pt(2), mid(2), mid(4)) + TriSurf(pt(2), mid(4), mid(1))
- + TriSurf(pt(3), mid(3), mid(4)) + TriSurf(pt(3), mid(4), mid(2));
- break;
- default :
- printf("Order %d quadrangle function space not implemented ", getPolynomialOrder()); print();
+const polynomialBasis* DI_Quad::getFunctionSpace(int o) const{
+ int order = (o == -1) ? getPolynomialOrder() : o;
+ switch (order) {
+ case 0: return polynomialBases::find(MSH_QUA_1);
+ case 1: return polynomialBases::find(MSH_QUA_4);
+ case 2: return polynomialBases::find(MSH_QUA_9);
+ case 3: return polynomialBases::find(MSH_QUA_16);
+ case 4: return polynomialBases::find(MSH_QUA_25);
+ case 5: return polynomialBases::find(MSH_QUA_36);
+ case 6: return polynomialBases::find(MSH_QUA_49);
+ case 7: return polynomialBases::find(MSH_QUA_64);
+ case 8: return polynomialBases::find(MSH_QUA_81);
+ case 9: return polynomialBases::find(MSH_QUA_100);
+ case 10: return polynomialBases::find(MSH_QUA_121);
+ default: Msg::Error("Order %d quadrangle function space not implemented", order);
}
}
-void DI_Quad::getShapeFunctions (double u, double v, double w, double s[], int ord) const {
- int order = (ord == -1) ? getPolynomialOrder() : ord;
- switch (order) {
- case 1 :
- s[0] = (1. - u) * (1. - v) / 4.;
- s[1] = (1. + u) * (1. - v) / 4.;
- s[2] = (1. + u) * (1. + v) / 4.;
- s[3] = (1. - u) * (1. + v) / 4.;
- break;
- case 2 :
- s[0] = u * v * (1. - u) * (1. - v) / 4.;
- s[1] = -u * v * (1. + u) * (1. - v) / 4.;
- s[2] = u * v * (1. + u) * (1. + v) / 4.;
- s[3] = -u * v * (1. - u) * (1. + v) / 4.;
- s[4] = -v * (1. - u) * (1. + u) * (1. - v) / 2.;
- s[5] = u * (1. + u) * (1. - v) * (1. + v) / 2.;
- s[6] = v * (1. - u) * (1. + u) * (1. + v) / 2.;
- s[7] = -u * (1. - u) * (1. - v) * (1. + v) / 2.;
- s[8] = (1. - u) * (1. + u) * (1. - v) * (1. + v);
- break;
- default :
- printf("Order %d quadrangle function space not implemented ", order); print();
- }
-}
-void DI_Quad::getGradShapeFunctions (const double u, const double v, const double w, double s[][3], int ord) const {
- int order = (ord == -1) ? getPolynomialOrder() : ord;
- switch (order) {
- case 1 :
- s[0][0] = -0.25 * (1. - v); s[0][1] = -0.25 * (1. - u); s[0][2] = 0.;
- s[1][0] = 0.25 * (1. - v); s[1][1] = -0.25 * (1. + u); s[1][2] = 0.;
- s[2][0] = 0.25 * (1. + v); s[2][1] = 0.25 * (1. + u); s[2][2] = 0.;
- s[3][0] = -0.25 * (1. + v); s[3][1] = 0.25 * (1. - u); s[3][2] = 0.;
- break;
- case 2 :
- s[0][0] = v * (1. - v) * (1. - 2. * u) / 4.;
- s[0][1] = u * (1. - u) * (1. - 2.* v ) / 4.;
- s[0][2] = 0.;
- s[1][0] = -v * (1. - v) * (1. + 2. * u) / 4.;
- s[1][1] = -u * (1. + u) * (1. - 2. * v) / 4.;
- s[1][2] = 0.;
- s[2][0] = v * (1. + v) * (1. + 2. * u) / 4.;
- s[2][1] = u * (1. + u) * (1. + 2. * v) / 4.;
- s[2][2] = 0.;
- s[3][0] = -v * (1. + v) * (1. - 2. * u) / 4.;
- s[3][1] = -u * (1. - u) * (1. + 2. * v) / 4.;
- s[3][2] = 0.;
- s[4][0] = v * (1. - v) * u;
- s[4][1] = -(1. - u * u) * (1. - 2. * v) / 2.;
- s[4][2] = 0.;
- s[5][0] = (1. - v * v) * (1. + 2. * u) / 2.;
- s[5][1] = -u * (1. + u) * v;
- s[5][2] = 0.;
- s[6][0] = -v * (1. + v) * u;
- s[6][1] = (1. - u * u) * (1. + 2. * v) / 2.;
- s[6][2] = 0.;
- s[7][0] = -(1. - v * v) * (1. - 2. * u) / 2.;
- s[7][1] = u * (1. - u) * v;
- s[7][2] = 0.;
- s[8][0] = -(1. - v * v) * u * 2.;
- s[8][1] = -(1. - u * u) * v * 2.;
- s[8][2] = 0.;
- break;
- default :
- printf("Order %d quadrangle function space not implemented ", order); print();
- }
+
+void DI_Quad::computeIntegral() {
+ integral_ = TriSurf(pt(0), pt(1), pt(2)) + TriSurf(pt(0), pt(2), pt(3));
+ // switch (getPolynomialOrder()) {
+ // case 1 :
+ // integral_ = TriSurf(pt(0), pt(1), pt(2)) + TriSurf(pt(0), pt(2), pt(3));
+ // break;
+ // case 2 :
+ // integral_ = TriSurf(pt(0), mid(0), mid(4)) + TriSurf(pt(0), mid(4), mid(3))
+ // + TriSurf(pt(1), mid(1), mid(4)) + TriSurf(pt(1), mid(4), mid(0))
+ // + TriSurf(pt(2), mid(2), mid(4)) + TriSurf(pt(2), mid(4), mid(1))
+ // + TriSurf(pt(3), mid(3), mid(4)) + TriSurf(pt(3), mid(4), mid(2));
+ // break;
+ // default :
+ // printf("Order %d quadrangle function space not implemented ", getPolynomialOrder()); print();
+ // }
}
// DI_Tetra methods -------------------------------------------------------------------------------
+const polynomialBasis* DI_Tetra::getFunctionSpace(int o) const{
+ int order = (o == -1) ? getPolynomialOrder() : o;
+ switch (order) {
+ case 0: return polynomialBases::find(MSH_TET_1);
+ case 1: return polynomialBases::find(MSH_TET_4);
+ case 2: return polynomialBases::find(MSH_TET_10);
+ case 3: return polynomialBases::find(MSH_TET_20);
+ case 4: return polynomialBases::find(MSH_TET_35);
+ case 5: return polynomialBases::find(MSH_TET_56);
+ case 6: return polynomialBases::find(MSH_TET_84);
+ case 7: return polynomialBases::find(MSH_TET_120);
+ case 8: return polynomialBases::find(MSH_TET_165);
+ case 9: return polynomialBases::find(MSH_TET_220);
+ case 10: return polynomialBases::find(MSH_TET_286);
+ default: Msg::Error("Order %d tetrahedron function space not implemented", order);
+ }
+}
+
void DI_Tetra::computeIntegral() {
integral_ = TetraVol(pt(0), pt(1), pt(2), pt(3));
}
-void DI_Tetra::getShapeFunctions (double u, double v, double w, double s[], int ord) const {
- int order = (ord == -1) ? getPolynomialOrder() : ord;
- switch (order) {
- case 1 :
- s[0] = 1. - u - v - w;
- s[1] = u;
- s[2] = v;
- s[3] = w;
- break;
- case 2 :
- s[0] = (1. - u - v - w) * (1. - 2. * u - 2. * v - 2. * w);
- s[1] = u * (2. * u - 1.);
- s[2] = v * (2. * v - 1.);
- s[3] = w * (2. * w - 1.);
- s[4] = 4. * u * (1. - u - v - w);
- s[5] = 4. * v * (1. - u - v - w);
- s[6] = 4. * w * (1. - u - v - w);
- s[7] = 4. * u * v;
- s[8] = 4. * v * w;
- s[9] = 4. * u * w;
- break;
- default :
- printf("Order %d tetrahedron function space not implemented ", order); print();
- }
-}
-void DI_Tetra::getGradShapeFunctions (const double u, const double v, const double w, double s[][3], int ord) const {
- int order = (ord == -1) ? getPolynomialOrder() : ord;
- switch (order) {
- case 1 :
- s[0][0] = -1.; s[0][1] = -1.; s[0][2] = -1.;
- s[1][0] = 1.; s[1][1] = 0.; s[1][2] = 0.;
- s[2][0] = 0.; s[2][1] = 1.; s[2][2] = 0.;
- s[3][0] = 0.; s[3][1] = 0.; s[3][2] = 1.;
- break;
- case 2 :
- s[0][0] = -3. + 4. * u + 4. * v + 4. * w;
- s[0][1] = -3. + 4. * u + 4. * v + 4. * w;
- s[0][2] = -3. + 4. * u + 4. * v + 4. * w;
- s[1][0] = 4. * u - 1.;
- s[1][1] = 0.;
- s[1][2] = 0.;
- s[2][0] = 0.;
- s[2][1] = 4. * v - 1.;
- s[2][2] = 0.;
- s[3][0] = 0.;
- s[3][1] = 0.;
- s[3][2] = 4. * w - 1.;
- s[4][0] = 4. * (1. - 2. * u - v - w);
- s[4][1] = -4. * u;
- s[4][2] = -4. * u;
- s[5][0] = -4. * v;
- s[5][1] = 4. * (1. - u - 2. * v - w);
- s[5][2] = -4. * v;
- s[6][0] = -4. * w;
- s[6][1] = -4. * w;
- s[6][2] = 4. * (1. - u - v - 2. * w);
- s[7][0] = 4. * v;
- s[7][1] = 4. * u;
- s[7][2] = 0.;
- s[8][0] = 0.;
- s[8][1] = 4. * w;
- s[8][2] = 4. * v;
- s[9][0] = 4. * w;
- s[9][1] = 0.;
- s[9][2] = 4. * u;
- break;
- default :
- printf("Order %d tetrahedron function space not implemented ", order); print();
- }
-}
double DI_Tetra::quality() const {
return qualityTet(x(0), y(0), z(0), x(1), y(1), z(1), x(2), y(2), z(2), x(3), y(3), z(3));
}
+
// Hexahedron methods -----------------------------------------------------------------------------
+const polynomialBasis* DI_Hexa::getFunctionSpace(int o) const{
+ int order = (o == -1) ? getPolynomialOrder() : o;
+ 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);
+ }
+}
void DI_Hexa::computeIntegral() {
integral_ = TetraVol(pt(0), pt(1), pt(3), pt(4)) + TetraVol(pt(1), pt(4), pt(5), pt(7))
+ TetraVol(pt(1), pt(3), pt(4), pt(7)) + TetraVol(pt(2), pt(5), pt(6), pt(7))
+ TetraVol(pt(1), pt(2), pt(3), pt(7)) + TetraVol(pt(1), pt(5), pt(2), pt(7));
}
-void DI_Hexa::getShapeFunctions (double u, double v, double w, double s[], int ord) const {
- int order = (ord == -1) ? getPolynomialOrder() : ord;
- switch (order) {
- case 1 :
- s[0] = (1. - u) * (1. - v) * (1. - w) / 8.;
- s[1] = (1. + u) * (1. - v) * (1. - w) / 8.;
- s[2] = (1. + u) * (1. + v) * (1. - w) / 8.;
- s[3] = (1. - u) * (1. + v) * (1. - w) / 8.;
- s[4] = (1. - u) * (1. - v) * (1. + w) / 8.;
- s[5] = (1. + u) * (1. - v) * (1. + w) / 8.;
- s[6] = (1. + u) * (1. + v) * (1. + w) / 8.;
- s[7] = (1. - u) * (1. + v) * (1. + w) / 8.;
- break;
- case 2 :
- s[0] = -u * v * w * (1. - u) * (1. - v) * (1. - w) / 8.;
- s[1] = u * v * w * (1. + u) * (1. - v) * (1. - w) / 8.;
- s[2] = -u * v * w * (1. + u) * (1. + v) * (1. - w) / 8.;
- s[3] = u * v * w * (1. - u) * (1. + v) * (1. - w) / 8.;
- s[4] = u * v * w * (1. - u) * (1. - v) * (1. + w) / 8.;
- s[5] = -u * v * w * (1. + u) * (1. - v) * (1. + w) / 8.;
- s[6] = u * v * w * (1. + u) * (1. + v) * (1. + w) / 8.;
- s[7] = -u * v * w * (1. - u) * (1. + v) * (1. + w) / 8.;
- s[8] = v * w * (1. - u) * (1. + u) * (1. - v) * (1. - w) / 4.;
- s[9] = -u * w * (1. - v) * (1. + v) * (1. + u) * (1. - w) / 4.;
- s[10] = -v * w * (1. - u) * (1. + u) * (1. + v) * (1. - w) / 4.;
- s[11] = u * w * (1. - v) * (1. + v) * (1. - u) * (1. - w) / 4.;
- s[12] = u * v * (1. - w) * (1. + w) * (1. - u) * (1. - v) / 4.;
- s[13] = -u * v * (1. - w) * (1. + w) * (1. + u) * (1. - v) / 4.;
- s[14] = u * v * (1. - w) * (1. + w) * (1. + u) * (1. + v) / 4.;
- s[15] = -u * v * (1. - w) * (1. + w) * (1. - u) * (1. + v) / 4.;
- s[16] = -v * w * (1. - u) * (1. + u) * (1. - v) * (1. + w) / 4.;
- s[17] = u * w * (1. - v) * (1. + v) * (1. + u) * (1. + w) / 4.;
- s[18] = v * w * (1. - u) * (1. + u) * (1. + v) * (1. + w) / 4.;
- s[19] = -u * w * (1. - v) * (1. + v) * (1. - u) * (1. + w) / 4.;
- s[20] = -w * (1. - w) * (1. + u) * (1. - u) * (1. + v) * (1. - v) / 2.;
- s[21] = -v * (1. - v) * (1. + u) * (1. - u) * (1. + w) * (1. - w) / 2.;
- s[22] = u * (1. + u) * (1. + v) * (1. - v) * (1. + w) * (1. - w) / 2.;
- s[23] = v * (1. + v) * (1. + u) * (1. - u) * (1. + w) * (1. - w) / 2.;
- s[24] = -u * (1. - u) * (1. + v) * (1. - v) * (1. + w) * (1. - w) / 2.;
- s[25] = w * (1. + w) * (1. + u) * (1. - u) * (1. + v) * (1. - v) / 2.;
- s[26] = (1. + u) * (1. - u) * (1. + v) * (1. - v) * (1. + w) * (1. - w);
- break;
- default :
- printf("Order %d hexahedron function space not implemented ", order); print();
- }
-}
-void DI_Hexa::getGradShapeFunctions (const double u, const double v, const double w, double s[][3], int ord) const {
- int order = (ord == -1) ? getPolynomialOrder() : ord;
- switch (order) {
- case 1 :
- s[0][0] = -0.125 * (1. - v) * (1. - w);
- s[0][1] = -0.125 * (1. - u) * (1. - w);
- s[0][2] = -0.125 * (1. - u) * (1. - v);
- s[1][0] = 0.125 * (1. - v) * (1. - w);
- s[1][1] = -0.125 * (1. + u) * (1. - w);
- s[1][2] = -0.125 * (1. + u) * (1. - v);
- s[2][0] = 0.125 * (1. + v) * (1. - w);
- s[2][1] = 0.125 * (1. + u) * (1. - w);
- s[2][2] = -0.125 * (1. + u) * (1. + v);
- s[3][0] = -0.125 * (1. + v) * (1. - w);
- s[3][1] = 0.125 * (1. - u) * (1. - w);
- s[3][2] = -0.125 * (1. - u) * (1. + v);
- s[4][0] = -0.125 * (1. - v) * (1. + w);
- s[4][1] = -0.125 * (1. - u) * (1. + w);
- s[4][2] = 0.125 * (1. - u) * (1. - v);
- s[5][0] = 0.125 * (1. - v) * (1. + w);
- s[5][1] = -0.125 * (1. + u) * (1. + w);
- s[5][2] = 0.125 * (1. + u) * (1. - v);
- s[6][0] = 0.125 * (1. + v) * (1. + w);
- s[6][1] = 0.125 * (1. + u) * (1. + w);
- s[6][2] = 0.125 * (1. + u) * (1. + v);
- s[7][0] = -0.125 * (1. + v) * (1. + w);
- s[7][1] = 0.125 * (1. - u) * (1. + w);
- s[7][2] = 0.125 * (1. - u) * (1. + v);
- break;
- default :
- printf("Order %d hexahedron function space not implemented ", order); print();
- }
-}
// ----------------------------------------------------------------------------
// -------------------------- SPLITTING ---------------------------------------
diff --git a/contrib/DiscreteIntegration/Integration3D.h b/contrib/DiscreteIntegration/Integration3D.h
index 22d3a6dcdab813b4fdbad4d454c85465a133a479..0748847f77744a4eec7cfe42a83730bebae38fb0 100644
--- a/contrib/DiscreteIntegration/Integration3D.h
+++ b/contrib/DiscreteIntegration/Integration3D.h
@@ -7,6 +7,8 @@
#include <map>
#include <cmath>
#include "gmshLevelset.h"
+#include "PolynomialBasis.h"
+#include "GmshDefines.h"
// Element type
#define DI_LIN 1
@@ -43,6 +45,16 @@ class DI_Point
DI_Point (double x, double y, double z, const double ls) : x_(x), y_(y), z_(z) {addLs(ls);}
DI_Point (double x, double y, double z, const DI_Element *e,
const std::vector<gLevelset *> &RPNi) : x_(x), y_(y), z_(z) {computeLs(e, RPNi);}
+ virtual const polynomialBasis* getFunctionSpace(int o) const
+ { return polynomialBases::find(MSH_PNT); }
+ virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
+ {
+ s[0] = 1.;
+ }
+ virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o)
+ {
+ s[0][0] = s[0][1] = s[0][2] = 0.;
+ }
// assignment
DI_Point & operator=(const DI_Point & rhs);
// add a levelset value (adjusted to 0 if ls<ZERO_LS_TOL) into the vector Ls
@@ -230,6 +242,7 @@ class DI_Element
if(pts_) delete [] pts_;
if(mid_) delete [] mid_;
}
+ virtual const polynomialBasis* getFunctionSpace(int order=-1) const { return 0; }
// return type
virtual int type() const = 0;
// return the dimension of the element
@@ -324,16 +337,17 @@ class DI_Element
inline int sizeLs() const {return pts_[0].sizeLs();}
// return the interpolating nodal shape functions evaluated at point (u,v,w)
// in parametric coordinates (if order = -1, use the polynomial order of the element)
- virtual void getShapeFunctions (double u, double v, double w, double s[], int order = -1) const = 0;
+ virtual void getShapeFunctions (double u, double v, double w, double s[], int order = -1) const;
+ // return the gradient of the shape functions evaluated at point (u,v,w) in the reference element
+ virtual void getGradShapeFunctions (const double u, const double v, const double w,
+ double s[][3], int order = -1) const;
// compute the coordinates in the element from the local coordinates (x,y,z)
void evalC (const double x, const double y, const double z, double *ev, int order = -1) const;
// evaluate the levelset at the local coordinates
// with the ith levelset value in the vector Ls of the points
// if i=-1, use the last value in Ls
double evalLs (const double x, const double y, const double z, int iLs = -1, int order = -1) const;
- // return the gradient of the shape functions evaluated at point (u,v,w) in the reference element
- virtual void getGradShapeFunctions (const double u, const double v, const double w,
- double s[][3], int order = -1) const = 0;
+
// compute the determinant of the jacobian at the point (u,v,w) in the reference element
double detJ (const double u, const double v, const double w) const;
// compute DetJ at each point, return false if 2 values have opposite sign
@@ -399,6 +413,7 @@ class DI_Line : public DI_Element
}
}
inline int nbEdg() const {return 1;}
+ virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
void computeIntegral();
inline double refIntegral() const {return 2.;}
inline DI_Point* pt (int i) const {return i < 2 ? &pts_[i] : &mid_[i - 2];}
@@ -414,9 +429,6 @@ class DI_Line : public DI_Element
void midV (const int e, int *s, int &n) const {
s[0] = 0; s[1] = 1; n = 2;
}
- void getShapeFunctions (double u, double v, double w, double s[], int order = -1) const;
- void getGradShapeFunctions (const double u, const double v, const double w,
- double s[][3], int order = -1) const;
double detJ (const double &xP, const double &yP, const double &zP) const;
bool cut (std::vector<gLevelset *> &LsRPN, std::vector<DI_IntegrationPoint *> &ip,
DI_Point::Container &cp, const int polynomialOrderL,
@@ -477,6 +489,7 @@ class DI_Triangle : public DI_Element
}
}
inline int nbEdg() const {return 3;}
+ virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
void computeIntegral();
inline double refIntegral() const {return 0.5;}
inline DI_Point* pt (int i) const {return i < 3 ? &pts_[i] : &mid_[i - 3];}
@@ -499,9 +512,6 @@ class DI_Triangle : public DI_Element
default : n = 0; return;
}
}
- void getShapeFunctions (double u, double v, double w, double s[], int order = -1) const;
- void getGradShapeFunctions (const double u, const double v, const double w,
- double s[][3], int order = -1) const;
double detJ (const double &xP, const double &yP, const double &zP) const;
bool cut (std::vector<gLevelset *> &LsRPN, std::vector<DI_IntegrationPoint *> &ip,
std::vector<DI_IntegrationPoint *> &ipS, DI_Point::Container &cp,
@@ -573,6 +583,7 @@ class DI_Quad : public DI_Element
}
}
inline int nbEdg() const {return 4;}
+ virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
void computeIntegral();
inline double refIntegral() const {return 4.;}
inline DI_Point* pt (int i) const {return i < 4 ? &pts_[i] : &mid_[i - 4];}
@@ -597,9 +608,6 @@ class DI_Quad : public DI_Element
default : n = 0; return;
}
}
- void getShapeFunctions (double u, double v, double w, double s[], int order = -1) const;
- void getGradShapeFunctions (const double u, const double v, const double w,
- double s[][3], int order = -1) const;
double detJ (const double &xP, const double &yP, const double &zP) const;
bool cut (std::vector<gLevelset *> &LsRPN, std::vector<DI_IntegrationPoint *> &ip,
std::vector<DI_IntegrationPoint *> &ipS, DI_Point::Container &cp,
@@ -667,6 +675,7 @@ class DI_Tetra : public DI_Element
}
}
inline int nbEdg() const {return 6;}
+ virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
void computeIntegral();
inline double refIntegral() const {return 1. / 6.;}
inline DI_Point* pt (int i) const {return i < 4 ? &pts_[i] : &mid_[i - 4];}
@@ -692,9 +701,6 @@ class DI_Tetra : public DI_Element
default : n = 0; return;
}
}
- void getShapeFunctions (double u, double v, double w, double s[], int order = -1) const;
- void getGradShapeFunctions (const double u, const double v, const double w,
- double s[][3], int order = -1) const;
double detJ (const double &xP, const double &yP, const double &zP) const;
bool cut (std::vector<gLevelset *> &LsRPN, std::vector<DI_IntegrationPoint *> &ip,
std::vector<DI_IntegrationPoint *> &ipS, DI_Point::Container &cp,
@@ -781,6 +787,7 @@ class DI_Hexa : public DI_Element
}
}
inline int nbEdg() const {return 12;}
+ virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
void computeIntegral();
inline double refIntegral() const {return 8.;}
inline DI_Point* pt (int i) const {return i < 8 ? &pts_[i] : &mid_[i - 8];}
@@ -821,9 +828,6 @@ class DI_Hexa : public DI_Element
default : n = 0; return;
}
}
- void getShapeFunctions (double u, double v, double w, double s[], int order = -1) const;
- void getGradShapeFunctions (const double u, const double v, const double w,
- double s[][3], int order = -1) const;
double detJ (const double &xP, const double &yP, const double &zP) const;
bool cut (std::vector<gLevelset *> &LsRPN, std::vector<DI_IntegrationPoint *> &ip,
std::vector<DI_IntegrationPoint *> &ipS, DI_Point::Container &cp,