diff --git a/Mesh/yamakawa.cpp b/Mesh/yamakawa.cpp
index a04ccd22381657c923e98045ce1028d2f52b14ca..f2b91b81295d392b03ba42e96d5e66d9abd53c97 100644
--- a/Mesh/yamakawa.cpp
+++ b/Mesh/yamakawa.cpp
@@ -30,6 +30,7 @@
#include <cassert>
+
#ifdef EXTERNAL_MIN_WEIGHTED_SET_SEARCH
#include "mwis.hpp"
#endif
@@ -5211,7 +5212,7 @@ bool PostOp::valid(MPyramid *pyr) {
return true;
}
-
+// Why template?
template <class T>
void export_the_clique_graphviz_format(cliques_compatibility_graph<T> &cl,
int clique_number, string filename)
@@ -5219,7 +5220,7 @@ void export_the_clique_graphviz_format(cliques_compatibility_graph<T> &cl,
ofstream out(filename.c_str());
out << "Graph G {" << endl;
typename multimap<int, set<T> >::reverse_iterator it_all = cl.allQ.rbegin();
- // multimap<int,set<T> >::reverse_iterator it_allen = cl.allQ.rend();
+
for (int i = 0; i < clique_number; i++) {
it_all++;
}
@@ -5227,9 +5228,6 @@ void export_the_clique_graphviz_format(cliques_compatibility_graph<T> &cl,
typename set<T>::iterator ithex = it_all->second.begin();
typename set<T>::iterator ithexen = it_all->second.end();
- //typedef tr1::unordered_multimap<hash_key, T> graph_data;
- //typedef tr1::unordered_multimap<hash_key, pair<T, graph_data > > graph;
-
int counter = 1;
map<T, int> visited_hex;
multimap<int, int> done;
@@ -5239,7 +5237,6 @@ void export_the_clique_graphviz_format(cliques_compatibility_graph<T> &cl,
typename cliques_compatibility_graph<T>::graph_data::const_iterator itgraphdata;
for (; itgraph != cl.end_graph(); itgraph++) {
-
T firsthex = itgraph->second.first;
// if (!post_check_validation(firsthex)) continue;
@@ -5249,9 +5246,9 @@ void export_the_clique_graphviz_format(cliques_compatibility_graph<T> &cl,
num1 = counter;
visited_hex[firsthex] = counter++;
}
- else
+ else {
num1 = itfind->second;
-
+ }
itgraphdata = itgraph->second.second.begin();
for (; itgraphdata != itgraph->second.second.end(); itgraphdata++) {
T secondhex = itgraphdata->second;
@@ -5262,9 +5259,9 @@ void export_the_clique_graphviz_format(cliques_compatibility_graph<T> &cl,
num2 = counter;
visited_hex[secondhex] = counter++;
}
- else
+ else {
num2 = itfind->second;
-
+ }
// search if num1 - num2 has already been written...
bool found = false;
pair<multimap<int, int>::iterator, multimap<int, int>::iterator> range =
@@ -5275,16 +5272,12 @@ void export_the_clique_graphviz_format(cliques_compatibility_graph<T> &cl,
break;
}
}
-
if (!found) {
done.insert(make_pair(num1, num2));
done.insert(make_pair(num2, num1));
out << num1 << " -- " << num2 << " ;" << endl;
}
-
}
-
-
}
// export chosen hex with different color
for (; ithex != ithexen; ithex++) { // brutal post-check: random pickup of hexahedra in clique
@@ -5611,7 +5604,14 @@ void cliques_compatibility_graph<T>::store_clique(int n)
};
-
+// Recursive function to find all maximal cliques in the input subgraph.
+// Implements more or less the algorithm proposed in
+// "The worst-case complexity for generating all maximal cliques and computational experiments"
+// Tomita et al., 2006. doi 10.1016/j.tcs.2006.06.015
+//
+// Very very very inefficient.
+//
+// What is n? the index of the clique being explored??.
template<class T>
void cliques_compatibility_graph<T>::find_cliques(graph_data &subgraph, int n)
{
@@ -5632,27 +5632,28 @@ void cliques_compatibility_graph<T>::find_cliques(graph_data &subgraph, int n)
T u = NULL; // Let's face it, this is a pointer.
hash_key u_key = 0;
// Choose the pivot: the node of the subgraph that
- // has a maximum number of neighbors in the candidate set.
+ // has a maximum number of adjacent nodes in the set of candidates
// Maximize the number of black nodes.
choose_u(subgraph, u, u_key);
- // The non-neighbors of u in the subgraph
+ // The non-neighbors of u in the subgraph (but should be candidates)
graph_data white;
- // The neighbors of u in the subgraph
+ // The nodes adjacent to u in the subgraph (but should be candidates)
graph_data black;
split_set_BW(u, u_key, subgraph, white, black);
- // recursive loop
+ // Explore all candidate and branches from there
while (white.size() > 0) {
Q.insert(u);
if (n == 0) {
- // WTF is this supposed to do? What is n ???
total_nodes_number = std::max(total_nodes_number, (unsigned int)white.size());
position++;
}
find_cliques(black, n + 1);
- if (cancel_search) break;
+ if (cancel_search) {
+ break;
+ }
erase_entry(white, u, u_key);
erase_entry(subgraph, u, u_key);
@@ -5684,11 +5685,9 @@ void cliques_compatibility_graph<T>::erase_entry(graph_data &subgraph, T &u, has
}
}
-// On veut choisir le sommet u qui maximise la taille de
-// l'intersection entre subgraph et les voisins de u dans le graphe
+// Chooses the node which ahs the maximum number of neighbors in the subgraph.
// Does not implement exactly what's is the paper [Tomita et al. 2006] where u
-// is chosen to maximize the intersection of its neighbors in the CAND subset. JP.
-// WHY is it not done like that here?
+// is chosen to maximize the intersection of its neighbors in the CAND subset. No idea why. JP.
template<class T>
void cliques_compatibility_graph<T>::choose_u(const graph_data &subgraph, T &u, hash_key &u_key)
{
@@ -5909,7 +5908,6 @@ size_t PEQuadrangle::get_max_nb_vertices()const { return 4; };
// and add it to the potential hex.
// TODO -- WE NEED ASSERTIONS in release mode and a command line to print stuff
// TODO -- Implement clear validity checks - topological validity and geometrical validity
-
// WARNING this function seems to be a mess -- I expect it to do not work as expected. JP.
// Stuff are commented, the test to check if a potential exist seems does not work
// and is bypassed.
@@ -6339,31 +6337,18 @@ void Recombinator_Graph::execute(GRegion* gr) {
create_losses_graph(gr);
- for (unsigned int i = 0; i < gr->tetrahedra.size(); i++) {
- MTetrahedron *tet = gr->tetrahedra[i];
- if (tet_to_hex.find(tet) == tet_to_hex.end()) {
- std::cout << "Cannot recombine: " << tet->getNum() << "\n";
- }
- }
-
#ifndef EXTERNAL_MIN_WEIGHTED_SET_SEARCH
compute_hex_ranks();
-
// a criteria to stop when the whole domain is exclusively composed of hex
found_the_ultimate_max_clique = false;
clique_stop_criteria<Hex*> criteria(hex_to_tet, gr->tetrahedra.size());
-
-
- cliques_losses_graph<Hex*> cl(incompatibility_graph, hex_ranks, max_nb_cliques, hex_to_tet.size(), &criteria, export_the_clique_graphviz_format);
+ cliques_losses_graph<Hex*> cl(incompatibility_graph, hex_ranks, max_nb_cliques, hex_to_tet.size(),
+ &criteria, export_the_clique_graphviz_format);
cl.find_cliques();
- //cl.export_cliques();
-
found_the_ultimate_max_clique = cl.found_the_ultimate_max_clique;
-
-
int clique_number = 0;
if (graphfilename.empty()) graphfilename.assign("mygraph.dot");
//export_clique_graphviz_format(cl,1,"mygraph2.dot");
@@ -6500,13 +6485,15 @@ bool Recombinator_Graph::merge_hex(GRegion *gr, Hex *hex) {
return true;
}
+
+
// Why derive then ?? this is quite stupid
void Recombinator_Graph::merge(GRegion* gr) {
throw;
}
-void Recombinator_Graph::export_tets(set<MElement*> &tetset, Hex* hex, string s)
-{
+
+void Recombinator_Graph::export_tets(set<MElement*> &tetset, Hex* hex, string s) {
stringstream ss;
ss << s.c_str();
ss << "hexptr_";
@@ -7196,6 +7183,7 @@ void Recombinator_Graph::add_face_connectivity(MElement *tet, int i, int j, int
delete t;
}
+
void Recombinator_Graph::compute_hex_ranks_blossom()
{
create_faces_connectivity();
@@ -7208,7 +7196,7 @@ void Recombinator_Graph::compute_hex_ranks_blossom()
}
map<Hex*, vector<double> >::iterator itfind = hex_ranks.find(hex);
if (itfind == hex_ranks.end())
- hex_ranks.insert(make_pair(hex, vector<double>(1)));
+ hex_ranks.insert(make_pair(hex, vector<double>(3)));
hex_ranks[hex][0] = nb_faces_on_boundary;
hex_ranks[hex][1] = hex->get_quality();
@@ -7219,7 +7207,6 @@ void Recombinator_Graph::compute_hex_ranks_blossom()
hex->vertex_in_facet(f, 2), hex->vertex_in_facet(f, 3));
if (face_in_blossom_info) count_blossom++;
}
-
hex_ranks[hex][2] = count_blossom;
}
}
diff --git a/Mesh/yamakawa.h b/Mesh/yamakawa.h
index caa2b74cb134543a3965290c5cf5211d06420b25..2a4c24ad706b3eba3a362a2847b0b0979a5960f3 100644
--- a/Mesh/yamakawa.h
+++ b/Mesh/yamakawa.h
@@ -28,14 +28,14 @@ class Hex {
private:
double quality;
unsigned long long hash;
- std::vector<MVertex*> vertices_;
+ std::vector<MVertex*> vertices_;
private:
void set_hash(){
hash = 0.;
for (int i = 0; i < 8; ++i) {
hash += vertices_[i]->getNum();
- }
+ }
}
public:
@@ -47,24 +47,24 @@ public:
set_hash();
}
Hex(MVertex* a2, MVertex* b2, MVertex* c2, MVertex* d2, MVertex* e2, MVertex* f2, MVertex* g2, MVertex* h2) :
- quality(0.)
+ quality(0.)
{
- vertices_.push_back(a2);
- vertices_.push_back(b2);
- vertices_.push_back(c2);
- vertices_.push_back(d2);
- vertices_.push_back(e2);
- vertices_.push_back(f2);
- vertices_.push_back(g2);
+ vertices_.push_back(a2);
+ vertices_.push_back(b2);
+ vertices_.push_back(c2);
+ vertices_.push_back(d2);
+ vertices_.push_back(e2);
+ vertices_.push_back(f2);
+ vertices_.push_back(g2);
vertices_.push_back(h2);
set_hash();
}
-
+
~Hex() {};
double get_quality() const { return quality; }
void set_quality(double new_quality) { quality = new_quality; }
-
+
MVertex* getVertex(unsigned int i) const {
if (i < 8) {
return vertices_[i];
@@ -151,7 +151,7 @@ public:
c = c2;
compute_hash();
}
- bool same_vertices(const Facet& facet) const {
+ bool same_vertices(const Facet& facet) const {
bool c1 = (a == facet.get_a()) || (a == facet.get_b()) || (a == facet.get_c());
bool c2 = (b == facet.get_a()) || (b == facet.get_b()) || (b == facet.get_c());
bool c3 = (c == facet.get_a()) || (c == facet.get_b()) || (c == facet.get_c());
@@ -214,7 +214,7 @@ private:
public:
Tuple() : v1(NULL), v2(NULL), v3(NULL), element(NULL), gf(NULL), hash(0.) {}
Tuple(MVertex* a, MVertex* b, MVertex* c, MElement* element2, GFace* gf2)
- {
+ {
MVertex* tmp[3] = { a,b,c };
std::sort(tmp, tmp + 3);
v1 = tmp[0];
@@ -222,10 +222,10 @@ public:
v2 = tmp[2];
hash = a->getNum() + b->getNum() + c->getNum();
element = element2;
- gf = gf2;
+ gf = gf2;
}
Tuple(MVertex* a, MVertex* b, MVertex* c)
- :element(NULL), gf(NULL)
+ :element(NULL), gf(NULL)
{
MVertex* tmp[3] = { a,b,c };
std::sort(tmp, tmp + 3);
@@ -254,7 +254,7 @@ public:
};
-// Class in charge of answering connectivity requests on
+// Class in charge of answering connectivity requests on
// the input tetraedral mesh
class TetMeshConnectivity {
public:
@@ -313,7 +313,7 @@ public:
TetSet tmp;
tets_around_vertices(v0, v1, tmp);
const TetSet& elements2 = tets_around_vertex(v2);
- std::set_intersection(tmp.begin(), tmp.end(),
+ std::set_intersection(tmp.begin(), tmp.end(),
elements2.begin(), elements2.end(), std::inserter(result, result.end()));
}
@@ -410,7 +410,7 @@ protected:
// Check if the hex is valid and compatible with the
// previously built hexes and add it to the region
bool add_hex_to_region_if_valid(const Hex& hex);
-
+
// ---- Computation of potential hexes
virtual void clear_potential_hex_info() {
potential.clear();
@@ -441,7 +441,7 @@ protected:
bool faces_statuquo(const Hex&);
bool faces_statuquo(MVertex*, MVertex*, MVertex*, MVertex*);
- bool are_all_tets_free(const std::set<MElement*>& tets) const;
+ bool are_all_tets_free(const std::set<MElement*>& tets) const;
void mark_tets(const std::set<MElement*>& tets);
void clear_hash_tables() {
@@ -467,7 +467,7 @@ protected:
// ---- Functions that should not be part of the class
double scaled_jacobian(MVertex*, MVertex*, MVertex*, MVertex*);
double max_scaled_jacobian(MElement*, int&);
- double min_scaled_jacobian(Hex&);
+ double min_scaled_jacobian(Hex&);
void print_statistics();
protected:
@@ -488,7 +488,7 @@ protected:
std::multiset<Diagonal> hash_tableC;
std::multiset<Tuple> tuples;
- std::set<MElement*> triangles;
+ std::set<MElement*> triangles;
};
@@ -531,11 +531,11 @@ public:
size_t get_max_nb_vertices() const;
};
-// Why are the following classes templates???
-// Used with only one class only, implemented in the cpp file (bad idea) and does not seem robust. JP.
+// Why are the following classes templates???
+// Used with only with Hex*, implemented in the cpp file (bad idea) and does not seem robust. JP.
// Very very complicated way to answer one question
-// Do we have all the tets of the input mesh for a given combination of hex?
+// Do we have all the tets of the input mesh for a given combination of hex?
// Slivers are tricky since they can be in shared by 2 hex.
template<class T>
class clique_stop_criteria {
@@ -551,9 +551,10 @@ private:
const unsigned int total_number_tet;
};
-// Why is this class template?
-// This complicates everything and is useless as the class
+// Why is this class template?
+// This complicates everything and is useless as the class
// cannot be used outside the .cpp file where it is implemented.
+// TODO - Rewrite this without multimaps or unnecessary abstractions.
template<class T>
class cliques_compatibility_graph {
public:
@@ -566,9 +567,9 @@ public:
typedef void(*ptrfunction_export)(cliques_compatibility_graph<T>&, int, string);
cliques_compatibility_graph(
- graph &_g,
+ graph &_g,
const map<T, std::vector<double> > &_hex_ranks,
- unsigned int _max_nb_cliques,
+ unsigned int _max_nb_cliques,
unsigned int _nb_hex_potentiels,
clique_stop_criteria<T> *csc,
ptrfunction_export fct);
@@ -581,10 +582,10 @@ public:
public:
bool found_the_ultimate_max_clique;
- // The stored maximal cliques.
- // The maximum number of stored cliques can be limited with the
+ // The stored maximal cliques.
+ // The maximum number of stored cliques can be limited with the
// max_nb_of_stored_cliques attribute
- // Cliques are ordered by their size (number of nodes in the clique)
+ // Cliques are ordered by size (number of nodes in the clique)
multimap<int, set<T> > allQ;
protected:
@@ -593,10 +594,9 @@ protected:
void split_set_BW(const T &u, const hash_key &u_key, const graph_data &subgraph, graph_data &white, graph_data &black);
void fill_black_set(const T &u, const hash_key &u_key, const graph_data &subgraph, graph_data &black);
void choose_u(const graph_data &subgraph, T &u, hash_key &u_key);
- // the maximum score (int) will be chosen...
double function_to_maximize_for_u(const T &u, const hash_key &u_key, const graph_data &subgraph);
void store_clique(int n);
- // returns true if two nodes are connected in the compatibility graph
+ // Returns true if two nodes are connected in the graph
virtual bool compatibility(const T &u, const hash_key &u_key, const T &v, const hash_key &v_key);
ptrfunction_export export_clique_graph;
@@ -612,12 +612,13 @@ protected:
unsigned int max_nb_of_stored_cliques;// to reduce memory footprint (set to zero if no limit)
clique_stop_criteria<T>* criteria;
bool cancel_search;
+ // Not used in anyway
const map<T, std::vector<double> > &hex_ranks;
graph &G;
graph_data_no_hash Q;// the current clique
};
-
+// Non necessary derivation
template<class T>
class cliques_losses_graph : public cliques_compatibility_graph<T> {
public:
@@ -627,21 +628,18 @@ public:
typedef void(*ptrfunction_export)(cliques_compatibility_graph<T>&, int, string);
cliques_losses_graph(
- graph &_g,
- const map<T, std::vector<double> > &_hex_ranks,
+ graph &_g,
+ const map<T, std::vector<double> > &_hex_ranks,
unsigned int _max_nb_cliques,
- unsigned int _nb_hex_potentiels,
- clique_stop_criteria<T> *csc,
+ unsigned int _nb_hex_potentiels,
+ clique_stop_criteria<T> *csc,
ptrfunction_export fct);
virtual ~cliques_losses_graph();
protected:
- // Returns true if the two nodes are compatible
+ // Returns true if the two nodes are compatible
// (connected in compatiblity graph - not connected in losses graph)
virtual bool compatibility(const T &u, const hash_key &u_key, const T &v, const hash_key &v_key);
- // AAAAaaaaahhhhhh exact same type and name for an attribute the
- // class from which this one derives.......
- // graph &G;
};
@@ -661,11 +659,10 @@ public:
typedef multimap<hash_key, Hex*> graph_data;
typedef multimap<hash_key, pair<Hex*, graph_data > > graph;
-
protected:
bool debug;
bool debug_graph;
-
+
int max_nb_cliques;
string graphfilename;
@@ -687,9 +684,9 @@ protected:
std::multimap<unsigned long long, Hex*> created_potential_hex;
- std::multimap<double, Hex*> degree;// degree = the final ranking of hexahedra
- std::multimap<int, Hex*> idegree;// idegree = number of connected hex in indirect neighbors graph
- std::multimap<int, Hex*> ndegree;// ndegree = number of direct neighbors !!! not chosen yet !!!
+ std::multimap<double, Hex*> degree; // degree = the final ranking of hexahedra
+ std::multimap<int, Hex*> idegree; // idegree = number of connected hex in indirect neighbors graph
+ std::multimap<int, Hex*> ndegree; // ndegree = number of direct neighbors !!! not chosen yet !!!
std::map<Hex*, int> reverse_idegree;
std::map<Hex*, int> reverse_ndegree;
// each tet has at least one neighbor, at most four. For all not chosen hex, check this data to find how many direct neighbors...
@@ -716,7 +713,9 @@ protected:
void fill_edges_table(const std::vector<MVertex*>&, Hex *hex);
void add_face(MVertex *a,MVertex* b,MVertex *c,Hex *hex);
void add_face(MVertex *a,MVertex* b,MVertex *c,std::multimap<unsigned long long, pair<PETriangle*,int> > &f);
-
+
+ // All the blossom related stuff is out of date - or not working
+ // Cannot be called. To remove?
bool find_face_in_blossom_info(MVertex *a, MVertex *b, MVertex *c, MVertex *d);
void compute_hex_ranks_blossom();
PETriangle* get_triangle(MVertex*a, MVertex* b, MVertex *c);
@@ -726,7 +725,7 @@ protected:
linemap::const_iterator find_the_line(PELine *t, const linemap &list);
std::multimap<unsigned long long, pair<PETriangle*,int> >::iterator
find_the_triangle(PETriangle *t, std::multimap<unsigned long long, pair<PETriangle*, int> > &list);
- std::multimap<unsigned long long, Hex* >::const_iterator
+ std::multimap<unsigned long long, Hex* >::const_iterator
find_the_created_potential_hex(Hex *t, const std::multimap<unsigned long long, Hex*> &list);
@@ -788,10 +787,12 @@ protected:
build_hash_tableC(hex);
}
- // Throw an assertion
+ // Throw an assertion
void merge(GRegion*);
// ------- exports --------
+ // ---- seems that it won't export nothing since the
+ // ---- data structures from which info is read seem to never be filled
void export_tets(set<MElement*> &tetset, Hex* hex, string s);
void export_single_hex_all(Hex* hex,string s);
void export_single_hex(Hex* hex,string s);
@@ -810,7 +811,7 @@ public:
virtual void buildGraphOnly(unsigned int max_nb_cliques, string filename=string());
virtual void buildGraphOnly(GRegion*, unsigned int max_nb_cliques, string filename=string());
virtual void execute_blossom(unsigned int max_nb_cliques, string filename=string());
- // What is this function supposed to do?
+ // What is this function supposed to do?
// Right now it throws at the first line. JP
virtual void execute_blossom(GRegion*, unsigned int max_nb_cliques, string filename=string());
virtual void createBlossomInfo();
@@ -971,7 +972,7 @@ public:
void trihedra(GRegion*);
void split_hexahedra(GRegion*);
void split_prisms(GRegion*);
- void split_pyramids(GRegion*);
+ void split_pyramids(GRegion*);
int nonConformDiag(MVertex* a,MVertex* b,MVertex* c,MVertex* d,GRegion* gr);
void pyramids1(MVertex*,MVertex*,MVertex*,MVertex*,GRegion*);
void pyramids2(MVertex*,MVertex*,MVertex*,MVertex*,GRegion*, bool allowNonConforming);
@@ -984,7 +985,7 @@ public:
//returns the geometrical validity of the pyramid
bool valid(MPyramid *pyr);
-
+
bool four(MElement*);
bool fourTrih(MElement*);
bool five(MElement*);