diff --git a/Geo/CMakeLists.txt b/Geo/CMakeLists.txt
index 4ed98e91e6e4d97f60a2de3fabfd815ca77a6c0e..b8763bb79e0f4b306e85930e005d68384d9e7ff0 100644
--- a/Geo/CMakeLists.txt
+++ b/Geo/CMakeLists.txt
@@ -34,6 +34,7 @@ set(SRC
SOrientedBoundingBox.cpp
GeomMeshMatcher.cpp
MVertex.cpp
+ MEdge.cpp
MFace.cpp
MElement.cpp MElementOctree.cpp
MLine.cpp MTriangle.cpp MQuadrangle.cpp MTetrahedron.cpp
diff --git a/Geo/MEdge.cpp b/Geo/MEdge.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..18a83bb5a7473caee99d9ea5b451ec045878f790
--- /dev/null
+++ b/Geo/MEdge.cpp
@@ -0,0 +1,63 @@
+// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#include "MEdge.h"
+#include "Numeric.h"
+
+bool MEdge::isInside(MVertex *v) const
+{
+ double tol = MVertexLessThanLexicographic::tolerance;
+ MVertex *v0 = _v[0];
+ MVertex *v1 = _v[1];
+ MVertexLessThanLexicographic lt;
+ if(lt(v0, v1)){
+ v0 = _v[1]; v1 = _v[0];
+ }
+ double x = v->x(), y = v->y(), z = v->z();
+ double x0 = v0->x(), y0 = v0->y(), z0 = v0->z();
+ double x1 = v1->x(), y1 = v1->y(), z1 = v1->z();
+ if(fabs(x - x0) < tol && fabs(y - y0) < tol && fabs(z - z0) < tol)
+ return true;
+ if(fabs(x - x1) < tol && fabs(y - y1) < tol && fabs(z - z1) < tol)
+ return true;
+ if(x < x0 - tol || x > x1 + tol ||
+ y < std::min(y0, y1) - tol || y > std::max(y0, y1) + tol ||
+ z < std::min(z0, z1) - tol || z > std::max(z0, z1) + tol)
+ ;return false;
+ if(fabs(x1 - x0) > tol){
+ double tx = (x - x0) / (x1 - x0);
+ if(fabs(y1 - y0) > tol){
+ double ty = (y - y0) / (y1 - y0);
+ if(fabs(z1 - z0) > tol){
+ double tz = (z - z0) / (z1 - z0);
+ if(fabs(tx - ty) > tol || fabs(tx - tz) > tol)
+ return false;
+ }
+ else{
+ if(fabs(tx - ty) > tol)
+ return false;
+ }
+ }
+ else{
+ if(fabs(z1 - z0) > tol){
+ double tz = (z - z0) / (z1 - z0);
+ if(fabs(tx - tz) > tol)
+ return false;
+ }
+ }
+ }
+ else{
+ if(fabs(y1 - y0) > tol){
+ double ty = (y - y0) / (y1 - y0);
+ if(fabs(z1 - z0) > tol){
+ double tz = (z - z0) / (z1 - z0);
+ if(fabs(ty - tz) > tol)
+ return false;
+ }
+ }
+ }
+ return true;
+}
+
diff --git a/Geo/MEdge.h b/Geo/MEdge.h
index c80a0cbba64e52bf1256da0192ec9efa27587dfa..700a788f9e079d2e5e046917e6450ec266dabde8 100644
--- a/Geo/MEdge.h
+++ b/Geo/MEdge.h
@@ -84,6 +84,7 @@ class MEdge {
t * _v[1]->y() + (1. - t) * _v[0]->y(),
t * _v[1]->z() + (1. - t) * _v[0]->z());
}
+ bool isInside(MVertex *v) const;
};
inline bool operator==(const MEdge &e1, const MEdge &e2)
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index 866d52238506717bbd26c04cdd4c09ba88d66342..6ef38440f58528c03a8041fd706cfa6cacb1fe7a 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -169,20 +169,51 @@ void MPolygon::_initVertices()
SVector3 ni = _parts[i]->getFace(0).normal();
if(dot(n, ni) < 0.) _parts[i]->revert();
}
-
- std::vector<MEdge> edg; // outer edges
+ // select only outer edges in vector edg
+ std::vector<MEdge> edg;
+ std::vector<MEdge> multiEdges;
for(int j = 0; j < _parts[0]->getNumEdges(); j++)
edg.push_back(_parts[0]->getEdge(j));
for(unsigned int i = 1; i < _parts.size(); i++) {
- for(int j = 0; j < 3; j++) {
+ for(int j = 0; j < _parts[i]->getNumEdges(); j++) {
+ bool found = false;
+ MEdge ed = _parts[i]->getEdge(j);
int k;
- for(k = edg.size() - 1; k >= 0; k--)
- if(_parts[i]->getEdge(j) == edg[k])
- break;
- if(k < 0)
- edg.push_back(_parts[i]->getEdge(j));
- else
- edg.erase(edg.begin() + k);
+ for(k = edg.size() - 1; k >= 0; k--){
+ if(ed == edg[k]){
+ edg.erase(edg.begin() + k);
+ found = true; break;
+ }
+ }
+ if(!found){
+ for(k = 0; k < multiEdges.size(); k++)
+ if(multiEdges[k].isInside(ed.getVertex(0)) &&
+ multiEdges[k].isInside(ed.getVertex(1))){
+ found = true; break;
+ }
+ }
+ if(!found){
+ for(k = edg.size() - 1; k >= 0; k--){
+ if(edg[k].isInside(ed.getVertex(0)) && edg[k].isInside(ed.getVertex(1))){
+ multiEdges.push_back(edg[k]);
+ edg.erase(edg.begin() + k);
+ found = true; break;
+ }
+ }
+ }
+ if(!found){
+ for(k = edg.size() - 1; k >= 0; k--){
+ if(ed.isInside(edg[k].getVertex(0)) && ed.isInside(edg[k].getVertex(1))){
+ edg.erase(edg.begin() + k);
+ int nbME = multiEdges.size();
+ if(nbME == 0 || multiEdges[nbME - 1] != ed)
+ multiEdges.push_back(ed);
+ found = true;
+ }
+ }
+ }
+ if(!found)
+ edg.push_back(ed);
}
}
diff --git a/contrib/DiscreteIntegration/Integration3D.cpp b/contrib/DiscreteIntegration/Integration3D.cpp
index 8219446b601dfe4e7f140807e64c0d1ef445542e..afa922c70eb17e293cf469a6163064818b050347 100644
--- a/contrib/DiscreteIntegration/Integration3D.cpp
+++ b/contrib/DiscreteIntegration/Integration3D.cpp
@@ -639,22 +639,23 @@ void DI_Element::setPolynomialOrder (int o) {
mid_ = NULL;
}
polOrder_ = o;
- switch (o) {
- case 1 :
- return;
- case 2 :
- mid_ = new DI_Point[nbMid()];
- for(int i = 0; i < nbMid(); i++) {
- std::vector<int> s(nbVert()); int n;
- midV(i, &s[0], n);
- double xc = 0, yc = 0, zc = 0;
- for(int j = 0; j < n; j++){
- xc += x(s[j]); yc += y(s[j]); zc += z(s[j]); }
- mid_[i] = DI_Point(xc/n, yc/n, zc/n);
- }
- return;
- default :
- printf("Order %d function space not implemented ", o); print();
+
+ if(o == 1) return;
+
+ const polynomialBasis *fs = getFunctionSpace(o);
+ if(!fs) Msg::Error("Function space not implemented for this type of element");
+
+ mid_ = new DI_Point[nbMid()];
+ int j = nbVert();
+ int dim = getDim();
+ double u, v, w;
+ double xyz[3];
+ for(int i = 0; i < nbMid(); i++, j++) {
+ u = fs->points(j,0);
+ v = (dim > 1) ? fs->points(j,1) : 0.;
+ w = (dim > 2) ? fs->points(j,2) : 0.;
+ evalC(u, v, w, xyz, 1);
+ mid_[i] = DI_Point(xyz[0], xyz[1], xyz[2]);
}
}
void DI_Element::setPolynomialOrder (int o, const DI_Element *e, const std::vector<gLevelset *> &RPNi) {
@@ -664,22 +665,23 @@ void DI_Element::setPolynomialOrder (int o, const DI_Element *e, const std::vect
mid_ = NULL;
}
polOrder_ = o;
- switch (o) {
- case 1 :
- return;
- case 2 :
- mid_ = new DI_Point[nbMid()];
- for(int i = 0; i < nbMid(); i++) {
- std::vector<int> s(nbVert()); int n;
- midV(i, &s[0], n);
- double xc = 0, yc = 0, zc = 0;
- for(int j = 0; j < n; j++){
- xc += x(s[j]); yc += y(s[j]); zc += z(s[j]); }
- mid_[i] = DI_Point(xc/n, yc/n, zc/n, e, RPNi);
- }
- return;
- default :
- printf("Order %d function space not implemented ", o); print();
+
+ if(o == 1) return;
+
+ const polynomialBasis *fs = getFunctionSpace(o);
+ if(!fs) Msg::Error("Function space not implemented for this type of element");
+
+ mid_ = new DI_Point[nbMid()];
+ int j = nbVert();
+ int dim = getDim();
+ double u, v, w;
+ double xyz[3];
+ for(int i = 0; i < nbMid(); i++, j++) {
+ u = fs->points(j,0);
+ v = (dim > 1) ? fs->points(j,1) : 0.;
+ w = (dim > 2) ? fs->points(j,2) : 0.;
+ evalC(u, v, w, xyz, 1);
+ mid_[i] = DI_Point(xyz[0], xyz[1], xyz[2], e, RPNi);
}
}
void DI_Element::addLs (const double *ls) {
@@ -1042,21 +1044,8 @@ bool DI_ElementLessThan::operator()(const DI_Element *e1, const DI_Element *e2)
// DI_Line methods --------------------------------------------------------------------------------
const polynomialBasis* DI_Line::getFunctionSpace(int o) const{
int order = (o == -1) ? getPolynomialOrder() : o;
- switch (order) {
- 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;
+ int tag = polynomialBasis::getTag(TYPE_LIN, order);
+ return polynomialBases::find(tag);
}
void DI_Line::computeIntegral() {
@@ -1076,20 +1065,8 @@ void DI_Line::computeIntegral() {
const polynomialBasis* DI_Triangle::getFunctionSpace(int o) const
{
int order = (o == -1) ? getPolynomialOrder() : o;
- switch (order) {
- 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);
- }
+ int tag = polynomialBasis::getTag(TYPE_TRI, order);
+ return polynomialBases::find(tag);
}
void DI_Triangle::computeIntegral() {
integral_ = TriSurf(pt(0), pt(1), pt(2));
@@ -1112,20 +1089,8 @@ double DI_Triangle::quality() const {
// DI_Quad methods --------------------------------------------------------------------------------
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);
- }
+ int tag = polynomialBasis::getTag(TYPE_QUA, order);
+ return polynomialBases::find(tag);
}
void DI_Quad::computeIntegral() {
@@ -1148,20 +1113,8 @@ void DI_Quad::computeIntegral() {
// 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);
- }
+ int tag = polynomialBasis::getTag(TYPE_TET, order);
+ return polynomialBases::find(tag);
}
void DI_Tetra::computeIntegral() {
@@ -1175,19 +1128,8 @@ double DI_Tetra::quality() const {
// 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);
- }
+ int tag = polynomialBasis::getTag(TYPE_HEX, order);
+ return polynomialBases::find(tag);
}
void DI_Hexa::computeIntegral() {
integral_ = TetraVol(pt(0), pt(1), pt(3), pt(4)) + TetraVol(pt(1), pt(4), pt(5), pt(7))