No commit message

No commit message

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]);
+	}
+	
+}
+

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();
     }