From cfc267eef990db45824ed8d1a3a05bf2de61b949 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Fri, 22 Nov 2013 19:18:33 +0000
Subject: [PATCH] oops again
---
Mesh/filterElements.cpp | 223 ++++++++++++++++++++++++++++++++++++++++
Mesh/filterElements.h | 23 +++++
2 files changed, 246 insertions(+)
create mode 100644 Mesh/filterElements.cpp
create mode 100644 Mesh/filterElements.h
diff --git a/Mesh/filterElements.cpp b/Mesh/filterElements.cpp
new file mode 100644
index 0000000000..f6f8473312
--- /dev/null
+++ b/Mesh/filterElements.cpp
@@ -0,0 +1,223 @@
+#include <algorithm>
+#include <vector>
+#include <set>
+#include "GmshDefines.h"
+#include "filterElements.h"
+#include "fullMatrix.h"
+#include "GaussIntegration.h"
+#include "MElement.h"
+#include "MTriangle.h"
+#include "MQuadrangle.h"
+#include "MPrism.h"
+#include "MHexahedron.h"
+#if defined(HAVE_RTREE)
+#include "rtree.h"
+void MElementBB(void *a, double *min, double *max);
+int MElementInEle(void *a, double *x);
+
+struct MElement_Wrapper
+{
+ bool _overlap;
+ MElement *_e;
+ std::vector<MElement*> _notOverlap;
+ MElement_Wrapper (MElement *e, std::vector<MElement*> notOverlap)
+ : _overlap(false), _e(e), _notOverlap(notOverlap)
+ {
+ std::sort(_notOverlap.begin(),_notOverlap.end());
+ }
+};
+
+static bool inBB (double *mi, double *ma, double *x){
+ if (x[0] < mi[0])return false;
+ if (x[1] < mi[1])return false;
+ if (x[2] < mi[2])return false;
+ if (x[0] > ma[0])return false;
+ if (x[1] > ma[1])return false;
+ if (x[2] > ma[2])return false;
+ return true;
+}
+
+bool rtree_callback(MElement *e1,void* pe2){
+ MElement_Wrapper *wrapper = static_cast<MElement_Wrapper*>(pe2);
+ MElement *e2 = wrapper->_e;
+
+ if (std::binary_search(wrapper->_notOverlap.begin(),wrapper->_notOverlap.end(),e1))return true;
+
+
+ for (int i=0;i<e1->getNumVertices();i++){
+ for (int j=0;j<e2->getNumVertices();j++){
+ if(e1->getVertex(i) == e2->getVertex(j))return true;
+ }
+ }
+
+ double min2[3],max2[3];
+ MElementBB(e2,min2,max2);
+
+ for (int i=0;i<e1->getNumVertices();i++){
+ SPoint3 xyz (e1->getVertex(i)->x(),e1->getVertex(i)->y(),e1->getVertex(i)->z());
+ if (inBB (min2,max2,xyz)){
+ if (MElementInEle(e2,xyz)){
+ wrapper->_overlap = true;
+ return false;
+ }
+ }
+ }
+
+ double min1[3],max1[3];
+ MElementBB(e1,min1,max1);
+ for (int i=0;i<e2->getNumVertices();i++){
+ SPoint3 xyz (e2->getVertex(i)->x(),e2->getVertex(i)->y(),e2->getVertex(i)->z());
+ if (inBB(min1,max1,xyz)){
+ if (MElementInEle(e1,xyz)){
+ wrapper->_overlap = true;
+ return false;
+ }
+ }
+ }
+
+ return false;
+
+ fullMatrix<double> pts1,pts2;
+ fullVector<double> weights;
+
+ gaussIntegration::get (e1->getType(),1,pts1,weights);
+ gaussIntegration::get (e2->getType(),1,pts2,weights);
+
+ for (int i=0;i<pts1.size1();i++){
+ double u = pts1(i,0);
+ double v = pts1(i,1);
+ double w = pts1(i,2);
+ SPoint3 xyz;
+ e1->pnt(u,v,w,xyz);
+ if (inBB(min2,max2,xyz)){
+ if (MElementInEle(e2,xyz)){
+ wrapper->_overlap = true;
+ return false;
+ }
+ }
+ }
+ for (int i=0;i<pts2.size1();i++){
+ double u = pts2(i,0);
+ double v = pts2(i,1);
+ double w = pts2(i,2);
+ SPoint3 xyz;
+ e2->pnt(u,v,w,xyz);
+ if (inBB(min1,max1,xyz)){
+ if (MElementInEle(e1,xyz)){
+ wrapper->_overlap = true;
+ return false;
+ }
+ }
+ }
+ return true;
+}
+
+struct Less_Partition : public std::binary_function<MElement*, MElement*, bool> {
+ bool operator()(const MElement *f1, const MElement *f2) const
+ {
+ return f1->getPartition() < f2->getPartition() ;
+ }
+};
+
+void filterColumns(std::vector<MElement*> &elem,
+ std::map<MElement*,std::vector<MElement*> > &_elemColumns)
+{
+ std::sort(elem.begin(),elem.end());
+ std::vector<MElement*> toKeep;
+ for (std::map<MElement*,std::vector<MElement*> >::iterator it = _elemColumns.begin();
+ it !=_elemColumns.end() ; ++it){
+ const std::vector<MElement*> &c = it->second;
+ unsigned int MAX = c.size() - 1;
+ for (unsigned int i=0;i<c.size(); i++){
+ if (!std::binary_search(elem.begin(),elem.end(),c[i])){
+ MAX = i - 1;
+ break;
+ }
+ }
+ // printf("MAX = %d c = %d\n",MAX,c.size());
+ for (unsigned int i=0;i<MAX;i++){
+ toKeep.push_back(c[i]);
+ }
+ for (unsigned int i=MAX;i<c.size();i++){
+ /// FIXME !!!
+ // delete c[i];
+ }
+ }
+ printf("%d --> %d\n",elem.size(),toKeep.size());
+ elem = toKeep;
+}
+
+
+void filterOverlappingElements (std::vector<MTriangle*> &blTris,
+ std::vector<MQuadrangle*>&blQuads,
+ std::map<MElement*,std::vector<MElement*> > &_elemColumns,
+ std::map<MElement*,MElement*> &_toFirst)
+{
+ std::vector<MElement*> vvv;
+ vvv.insert(vvv.begin(),blTris.begin(),blTris.end());
+ vvv.insert(vvv.begin(),blQuads.begin(),blQuads.end());
+ Less_Partition lp;
+ std::sort(vvv.begin(),vvv.end(), lp);
+ filterOverlappingElements (vvv,_elemColumns,_toFirst);
+ filterColumns (vvv,_elemColumns);
+ blTris.clear();
+ blQuads.clear();
+ for (unsigned int i=0;i<vvv.size();i++){
+ if (vvv[i]->getType() == TYPE_TRI)blTris.push_back((MTriangle*)vvv[i]);
+ else if (vvv[i]->getType() == TYPE_QUA)blQuads.push_back((MQuadrangle*)vvv[i]);
+ }
+}
+
+void filterOverlappingElements (std::vector<MPrism*> &blPrisms,
+ std::vector<MHexahedron*>&blHexes,
+ std::map<MElement*,std::vector<MElement*> > &_elemColumns,
+ std::map<MElement*,MElement*> &_toFirst)
+{
+ printf("filtering !!\n");
+ std::vector<MElement*> vvv;
+ vvv.insert(vvv.begin(),blPrisms.begin(),blPrisms.end());
+ vvv.insert(vvv.begin(),blHexes.begin(),blHexes.end());
+ Less_Partition lp;
+ std::sort(vvv.begin(),vvv.end(), lp);
+ filterOverlappingElements (vvv,_elemColumns,_toFirst);
+ filterColumns (vvv,_elemColumns);
+ blPrisms.clear();
+ blHexes.clear();
+ for (unsigned int i=0;i<vvv.size();i++){
+ if (vvv[i]->getType() == TYPE_PRI)blPrisms.push_back((MPrism*)vvv[i]);
+ else if (vvv[i]->getType() == TYPE_HEX)blHexes.push_back((MHexahedron*)vvv[i]);
+ }
+}
+
+
+void filterOverlappingElements (std::vector<MElement*> &els,
+ std::map<MElement*,std::vector<MElement*> > &_elemColumns,
+ std::map<MElement*,MElement*> &_toFirst )
+{
+ std::vector<MElement*> newEls;
+ RTree<MElement*,double,3,double> rtree;
+ for (unsigned int i=0;i<els.size();i++){
+ MElement *e = els[i];
+ double _min[3],_max[3];
+ MElementBB(e,_min,_max);
+ MElement *first = _toFirst[e];
+ MElement_Wrapper w(e,_elemColumns[first]);
+ rtree.Search(_min,_max,rtree_callback,&w);
+ if (w._overlap){
+ delete e;
+ }
+ else {
+ rtree.Insert(_min,_max,e);
+ newEls.push_back(e);
+ }
+ }
+ els = newEls;
+}
+#else
+void filterOverlappingElements (std::vector<MElement*> &els,
+ std::map<MElement*,std::vector<MElement*> > &_elemColumns,
+ std::map<MElement*,MElement*> &_toFirst )
+{
+ Msg::Error("full boundary layer capabilities are only available whilst compiling GMSH together with a RTREE library");
+}
+#endif
diff --git a/Mesh/filterElements.h b/Mesh/filterElements.h
new file mode 100644
index 0000000000..f046b37d63
--- /dev/null
+++ b/Mesh/filterElements.h
@@ -0,0 +1,23 @@
+#ifndef _FILTER_OVERLAPPING_ELEMENTS_
+#define _FILTER_OVERLAPPING_ELEMENTS_
+#include <map>
+#include <vector>
+class MElement;
+class MPrism;
+class MHexahedron;
+class MTriangle;
+class MQuadrangle;
+void filterOverlappingElements (std::vector<MElement*> &,
+ std::map<MElement*,std::vector<MElement*> > &_elemColumns,
+ std::map<MElement*,MElement*> &_toFirst);
+void filterOverlappingElements (std::vector<MPrism*> &blPrisms,
+ std::vector<MHexahedron*>&blHexes,
+ std::map<MElement*,std::vector<MElement*> > &_elemColumns,
+ std::map<MElement*,MElement*> &_toFirst);
+void filterOverlappingElements (std::vector<MTriangle*> &blTris,
+ std::vector<MQuadrangle*>&blQuads,
+ std::map<MElement*,std::vector<MElement*> > &_elemColumns,
+ std::map<MElement*,MElement*> &_toFirst);
+void filterColumns(std::vector<MElement*> &elem,
+ std::map<MElement*,std::vector<MElement*> > &_elemColumns);
+#endif
--
GitLab