diff --git a/Solver/filters.cpp b/Solver/filters.cpp
index d48467dc84152d5b427f52ae82e3645bf22b8a7e..9869181b16cf3e12a6f9a7c63d636035a0b7cbfe 100644
--- a/Solver/filters.cpp
+++ b/Solver/filters.cpp
@@ -9,8 +9,259 @@
//
#include "filters.h"
+#include "MTriangle.h"
+#include "MQuadrangle.h"
+#include "MLine.h"
+#include "MTetrahedron.h"
+#include "MHexahedron.h"
+#include "MElementCut.h"
+#include "MElement.h"
void FilterLevelSetForLagMultSpace::SortNodes (void)
{
- ;
-}
\ No newline at end of file
+
+ std::set< EdgeType > Se; // edges
+ std::set< NodeType > Sn; // level set nodes kept
+ std::map < NodeType, int > NodesScore;
+ std::map < NodeType , EdgeType > NodeToEdgeMap;
+ std::set< NodeType >::iterator itn;
+
+ // initialization
+ fillNodeToEdgeMap (NodeToEdgeMap,Se,Sn);
+
+ //-- decimation algorithm --
+
+
+ while (!Sn.empty())
+ {
+ // compute score of Sn
+ ComputeScore(NodeToEdgeMap,Se,Sn,NodesScore);
+ // take the lowest score
+ NodeType v = getNodeWithLowestScore(NodesScore);
+ // if NOT marked as looser then it s a winner
+ if (_looser_nodes.find(v) != _looser_nodes.end()) _winner_nodes.insert(v);
+ // kill all connected edges and nodes (then mark them as loosers)
+ killConnectedEdges (NodeToEdgeMap,Se,Sn,v);
+ }
+
+
+}
+
+void FilterLevelSetForLagMultSpace::fillNodeToEdgeMap(std::map < NodeType , EdgeType > & NodeToEdgeMap, std::set< EdgeType > & Se , std::set< NodeType > & Sn)
+{
+ std::set<MElement*>::const_iterator it = _LevelSetElements->begin();
+ for (; it != _LevelSetElements->end(); it++)
+ {
+ MElement *e = *it;
+ for (int i = 0 ; i < e->getNumVertices() ; i ++) // warning getnumvertices polygon : vertices + inner ?
+ {
+ NodeType v = e->getVertex(i);
+ EdgeType edge = findEdge(v,e);
+ Sn.insert(v);
+ Se.insert(edge);
+ NodeToEdgeMap[v] = edge;
+ }
+ }
+}
+
+std::pair <MVertex * , MVertex * > FilterLevelSetForLagMultSpace::findEdge(NodeType v, MElement * e)
+{
+ MElement *ep;
+ if (e->getParent()) ep = e->getParent();
+ switch (ep->getType())
+ {
+ case TYPE_TRI :
+ {
+ EdgeType edge;
+ for (int i = 0 ; i < 3 ; i ++) // for all edges
+ {
+ int n1,n2;
+ n1 = i;
+ n2 = (i+1)%3;
+ // if it s an element vertex
+ if (v->getNum() == ep->getVertex(n1)->getNum())
+ {
+ edge = EdgeType (ep->getVertex(n1),ep->getVertex(n1));
+ return edge;
+ }
+ if (v->getNum() == ep->getVertex(n2)->getNum())
+ {
+ edge = EdgeType (ep->getVertex(n2),ep->getVertex(n2));
+ return edge;
+ }
+ // else if it s not an element vertex
+ edge = EdgeType(ep->getVertex(n1),ep->getVertex(n2));
+ if (NodeBelongToEdge(v,edge)) return edge;
+ }
+ }
+ case TYPE_QUA :
+ {
+ EdgeType edge;
+ for (int i = 0 ; i < 4 ; i ++) // for all edges
+ {
+ int n1,n2;
+ n1 = i;
+ n2 = (i+1)%4;
+ // if it s an element vertex
+ if (v->getNum() == ep->getVertex(n1)->getNum())
+ {
+ edge = EdgeType (ep->getVertex(n1),ep->getVertex(n1));
+ return edge;
+ }
+ if (v->getNum() == ep->getVertex(n2)->getNum())
+ {
+ edge = EdgeType (ep->getVertex(n2),ep->getVertex(n2));
+ return edge;
+ }
+ // else if it s not an element vertex
+ edge = EdgeType(ep->getVertex(n1),ep->getVertex(n2));
+ if (NodeBelongToEdge(v,edge)) return edge;
+ }
+ }
+ case TYPE_TET :
+ {
+ int tab[6][2] = {{0, 1}, {0, 3}, {0, 2}, {1, 2}, {2, 3}, {3, 1}}; // edges
+ EdgeType edge;
+ for (int i = 0 ; i < 6 ; i ++) // for all edges
+ {
+ int n1,n2;
+ n1 = tab[i][0];
+ n2 = tab[i][1];
+ // if it s an element vertex
+ if (v->getNum() == ep->getVertex(n1)->getNum())
+ {
+ edge = EdgeType (ep->getVertex(n1),ep->getVertex(n1));
+ return edge;
+ }
+ if (v->getNum() == ep->getVertex(n2)->getNum())
+ {
+ edge = EdgeType (ep->getVertex(n2),ep->getVertex(n2));
+ return edge;
+ }
+ // else if it s not an element vertex
+ edge = EdgeType(ep->getVertex(n1),ep->getVertex(n2));
+ if (NodeBelongToEdge(v,edge)) return edge;
+ }
+ }
+ case TYPE_HEX :
+ {
+ int tab[12][2] = {{0, 1}, {0, 4}, {0, 3}, {1, 2}, {1, 5},
+ {2, 3},{2,6},{3, 7}, {4, 7}, {4, 5},{5, 6},{6, 7}};// edges
+ EdgeType edge;
+ for (int i = 0 ; i < 12 ; i ++) // for all edges
+ {
+ int n1,n2;
+ n1 = tab[i][0];
+ n2 = tab[i][1];
+ // if it s an element vertex
+ if (v->getNum() == ep->getVertex(n1)->getNum())
+ {
+ edge = EdgeType (ep->getVertex(n1),ep->getVertex(n1));
+ return edge;
+ }
+ if (v->getNum() == ep->getVertex(n2)->getNum())
+ {
+ edge = EdgeType (ep->getVertex(n2),ep->getVertex(n2));
+ return edge;
+ }
+ // else if it s not an element vertex
+ edge = EdgeType(ep->getVertex(n1),ep->getVertex(n2));
+ if (NodeBelongToEdge(v,edge)) return edge;
+ }
+ }
+ default : std::cout << "findEdge warning..." << std::endl;
+ }
+}
+
+
+bool FilterLevelSetForLagMultSpace::NodeBelongToEdge(NodeType v, EdgeType edge)
+{
+
+ double eps = 1e-10;
+
+ double distn, edgelength;
+
+ edgelength = edge.first->distance(edge.second);
+ distn = v->distance(edge.first) + v->distance(edge.second);
+
+ if (std::abs(edgelength - distn) < eps) return true;
+ else return false;
+
+}
+
+
+void FilterLevelSetForLagMultSpace::ComputeScore(
+std::map < NodeType , EdgeType > & NodeToEdgeMap, std::set< EdgeType > & Se, std::set< NodeType > & Sn, std::map < NodeType , int > & NodesScore)
+{
+ NodesScore.clear();
+ std::set < NodeType >::iterator it = Sn.begin();
+ for (;it!=Sn.end();it++)
+ {
+ int score;
+ EdgeType edge;
+ edge = NodeToEdgeMap[(*it)];
+ if (edge.first->getNum() == edge.second->getNum()) score = 0; // regular mesh nodes score
+ else score = ComputeIncidentEdges(Se,edge.first) + ComputeIncidentEdges(Se,edge.second);
+ NodesScore[(*it)] = score;
+ }
+}
+
+
+int FilterLevelSetForLagMultSpace::ComputeIncidentEdges(std::set< EdgeType > & Se, NodeType & v)
+{
+ int n = 0;
+ std::set < EdgeType >::iterator it = Se.begin();
+ for (;it!=Se.end();it++)
+ {
+ if ((*it).first == v | (*it).second == v) n++;
+ }
+ return n;
+}
+
+
+MVertex * FilterLevelSetForLagMultSpace::getNodeWithLowestScore(std::map < NodeType , int > & NodesScore)
+{
+ NodeType v;
+ std::map < NodeType , int >::iterator it = NodesScore.begin();
+ v = (*it).first;
+ for (;it!=NodesScore.end() ; it ++)
+ {
+ if ((*it).second < NodesScore[v]) v = (*it).first ;
+ }
+ return v;
+}
+
+void FilterLevelSetForLagMultSpace::killConnectedEdges (std::map < NodeType , EdgeType > & NodeToEdgeMap, std::set< EdgeType > & Se , std::set< NodeType > & Sn, NodeType v)
+{
+
+ std::vector < EdgeType > EdgesToErase;
+ std::vector < NodeType > NodesToErase;
+
+ NodeType v1 = NodeToEdgeMap[v].first;
+ NodeType v2 = NodeToEdgeMap[v].second;
+
+ std::set < EdgeType >::iterator ite = Se.begin();
+ for (;ite!=Se.end();ite++)
+ {
+ if ((*ite).first == v1 | (*ite).second == v1) EdgesToErase.push_back(*ite);
+ if ((*ite).first == v2 | (*ite).second == v2) EdgesToErase.push_back(*ite);
+ }
+ for (int i = 0 ; i < EdgesToErase.size() ; i ++)
+ {
+ Se.erase(EdgesToErase[i]);
+ }
+
+
+ std::set < NodeType >::iterator itn = Sn.begin();
+ for (; itn != Sn.end() ; itn ++)
+ {
+ if (Se.find(NodeToEdgeMap[(*itn)]) == Se.end()) NodesToErase.push_back(*itn);
+ }
+ for (int i = 0 ; i < NodesToErase.size() ; i ++)
+ {
+ if (NodesToErase[i]->getNum() != v->getNum()) _looser_nodes.insert(NodesToErase[i]);
+ Sn.erase(NodesToErase[i]);
+ }
+
+}
+
diff --git a/Solver/filters.h b/Solver/filters.h
index 9a5781143c3b0c038d58fd85b146e4390eb30ef4..74fa9615c40eb83765458e2b3eba5f941b8f40a3 100644
--- a/Solver/filters.h
+++ b/Solver/filters.h
@@ -15,6 +15,7 @@
#include "dofManager.h"
#include "GModel.h"
#include "groupOfElements.h"
+#include "DILevelset.h"
class FilterNodeEnriched
{
@@ -124,21 +125,45 @@ class FilterElementsCutByLevelSet
class FilterLevelSetForLagMultSpace
{
+ private :
+
+ typedef MVertex * NodeType;
+ typedef std::pair <NodeType , NodeType > EdgeType;
+ typedef std::pair <EdgeType , double > EdgeScoreType;
+
+
private :
- groupOfElements * _g;
+ groupOfElements * _LevelSetElements;
std::pair<int,int> _LevelSetEntity;
- std::set<int> _winner_nodes;
- std::set<int> _all_nodes;
+ std::set< NodeType > _winner_nodes;
+ std::set< NodeType > _looser_nodes;
+ gLevelset *_ls;
private :
+ // decimation algorithm result in _winner_nodes set
void SortNodes (void) ;
+ // initialisation of needed sets
+ void fillNodeToEdgeMap(std::map < NodeType , EdgeType > & NodeToEdgeMap, std::set< EdgeType > & Se , std::set< NodeType > & Sn) ;
+ // find edge in the element associate with node
+ EdgeType findEdge(NodeType v, MElement * e);
+ // verify if node belong to edge // ...normaly has to be done when cutting model ...
+ bool NodeBelongToEdge(NodeType v, EdgeType edge);
+ // compute score
+ void ComputeScore(std::map < NodeType , EdgeType > & NodeToEdgeMap, std::set< EdgeType > & Se, std::set< NodeType > & Sn, std::map < NodeType , int > & NodesScore);
+ // compute number of incident edges in Se to node v
+ int ComputeIncidentEdges(std::set< EdgeType > & Se, NodeType & v);
+ // get lowest (no need to sort just get the lowest)
+ NodeType getNodeWithLowestScore(std::map < NodeType , int > & NodesScore);
+ // kill connected edges
+ void killConnectedEdges (std::map < NodeType , EdgeType > & NodeToEdgeMap, std::set< EdgeType > & Se , std::set< NodeType > & Sn, NodeType v);
public :
- FilterLevelSetForLagMultSpace(std::pair<int,int> LevelSetEntity, groupOfElements * g) : _LevelSetEntity(LevelSetEntity), _g(g)
+ FilterLevelSetForLagMultSpace(std::pair<int,int> LevelSetEntity , gLevelset *ls) : _LevelSetEntity(LevelSetEntity), _ls(ls)
{
+ groupOfElements *_LevelSetElements = new groupOfElements (_LevelSetEntity.first, _LevelSetEntity.second);
SortNodes();
}