yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #01045
[svn] r1696 - trunk/extra/SpherePadder
Author: richefeu
Date: 2009-02-26 13:47:38 +0100 (Thu, 26 Feb 2009)
New Revision: 1696
Modified:
trunk/extra/SpherePadder/SpherePadder.cpp
trunk/extra/SpherePadder/SpherePadder.hpp
trunk/extra/SpherePadder/TetraMesh.cpp
trunk/extra/SpherePadder/main.cpp
Log:
Begin devel. of densification process
Modified: trunk/extra/SpherePadder/SpherePadder.cpp
===================================================================
--- trunk/extra/SpherePadder/SpherePadder.cpp 2009-02-26 08:35:54 UTC (rev 1695)
+++ trunk/extra/SpherePadder/SpherePadder.cpp 2009-02-26 12:47:38 UTC (rev 1696)
@@ -52,8 +52,8 @@
trace_functions = true;
meshIsPlugged = false;
probeIsDefined = false;
- ratio = 5.0; rmoy = 0.0;
- virtual_radius_factor = 5.0;
+ ratio = 2.0; rmoy = 0.0;
+ virtual_radius_factor = 10.0;
/* FIXME
pour le moment, l'utilisateur ne peut entre qu'un ratio.
@@ -171,9 +171,8 @@
place_at_tetra_vertexes ();
//save_granulo("tetra_vertexes_granulo.dat");
//detect_overlap();
- //place_virtual_spheres(); // TODO deplacer dans une fonction 'densify' par exemple
+ //place_virtual_spheres();
-
time_t stop_time = clock();
cerr << "Total number of spheres = " << sphere.size() << endl;
@@ -187,6 +186,37 @@
cerr << "Time used = " << time_used << " s" << endl;
}
+void SpherePadder::densify()
+{
+ //BEGIN_FUNCTION ("Densify");
+
+ //place_virtual_spheres();
+ for (unsigned int i = 0 ; i < sphere.size() ; i++)
+ {
+ triangulation.insert_node(sphere[i].x, sphere[i].y, sphere[i].z, i, false);
+ }
+
+ tetra_porosity P;
+ if (!tetra_porosities.empty()) tetra_porosities.clear();
+ // if triangulation is empty ...
+ triangulation.init_current_tetrahedron();
+ do
+ {
+ triangulation.current_tetrahedron_get_nodes(P.id1,P.id2,P.id3,P.id4);
+ P.volume = triangulation.current_tetrahedron_get_volume();
+ P.void_volume = P.volume - solid_volume_of_tetrahedron(sphere[P.id1], sphere[P.id2], sphere[P.id3], sphere[P.id4]);
+ tetra_porosities.push_back(P);
+ } while (triangulation.next_tetrahedron());
+ qsort(&(tetra_porosities[0]),tetra_porosities.size(),sizeof(tetra_porosity),compare_tetra_porosity);
+
+ for (unsigned int i = 0 ; i < tetra_porosities.size() ; i++)
+ {
+ cout << tetra_porosities[i].volume << "\t" << tetra_porosities[i].void_volume << endl;
+ }
+
+ //END_FUNCTION;
+}
+
void SpherePadder::tetra_pad() // EN TRAVAUX !!!!!!!!!!!!!!!!!!!!!!!!!
{
place_at_nodes();
@@ -277,16 +307,17 @@
// tmp (bricolage)
if (i < mesh->node.size())
+ {
for (unsigned int s = 0 ; s < mesh->segment.size() ; ++s)
- {
- if (mesh->segment[s].nodeId[0] == i)
- {
- if (mesh->segment[s].nodeId[1] < mesh->node.size())
- fmgpost << " <SPSPx antac=\"" << mesh->segment[s].nodeId[1] + 1 << "\"/>" << endl;
- }
- }
+ {
+ if (mesh->segment[s].nodeId[0] == i)
+ {
+ if (mesh->segment[s].nodeId[1] < mesh->node.size())
+ fmgpost << " <SPSPx antac=\"" << mesh->segment[s].nodeId[1] + 1 << "\"/>" << endl;
+ }
+ }
+ }
-
fmgpost << " </body>" << endl;
}
@@ -970,13 +1001,11 @@
k = sphere_virtual_radius;
sphereId = sphere.size();
-
for (unsigned int f = 0 ; f < mesh->face.size() ; ++f)
{
-
- if ( mesh->face[f].belongBoundary == true )
- {
-
+
+ if (mesh->face[f].belongBoundary)
+ {
current_tetra_id = mesh->face[f].tetraOwner[0];
current_tetra = mesh->tetraedre[current_tetra_id];
@@ -1139,7 +1168,7 @@
S.tetraOwner = current_tetra_id;
sphere.push_back(S);
partition.add(sphereId++,S.x,S.y,S.z);
-
+
/////////////////////////////////////////////
}
@@ -1415,7 +1444,7 @@
}
-
+// point : x,y,z,R
double SpherePadder::spherical_triangle (double point1[],double point2[],double point3[],double point4[])
{
@@ -1423,7 +1452,6 @@
if (rayon == 0.0) return 0.0;
double vect12[3], vect13[3], vect14[3];
- const double pi = 3.14159265;
vect12[0] = point2[0] - point1[0];
vect12[1] = point2[1] - point1[1];
@@ -1441,6 +1469,10 @@
double norme13 = sqrt( (vect13[0]*vect13[0]) + (vect13[1]*vect13[1]) + (vect13[2]*vect13[2]) );
double norme14 = sqrt( (vect14[0]*vect14[0]) + (vect14[1]*vect14[1]) + (vect14[2]*vect14[2]) );
+// if (norme12 == 0.0) cout << "sin(A) == 0.0" << endl;
+// if (norme13 == 0.0) cout << "sin(B) == 0.0" << endl;
+// if (norme14 == 0.0) cout << "sin(C) == 0.0" << endl;
+
double A = acos( (vect12[0]*vect13[0] + vect12[1]*vect13[1] + vect12[2]*vect13[2]) / (norme13 * norme12));
double B = acos( (vect12[0]*vect14[0] + vect12[1]*vect14[1] + vect12[2]*vect14[2]) / (norme14 * norme12));
double C = acos( (vect14[0]*vect13[0] + vect14[1]*vect13[1] + vect14[2]*vect13[2]) / (norme13 * norme14));
@@ -1449,14 +1481,16 @@
double b = acos( (cos(B) - cos(C) * cos(A)) / (sin(C) * sin(A)) );
double c = acos( (cos(C) - cos(A) * cos(B)) / (sin(A) * sin(B)) );
+// if (sin(A) == 0.0) cout << "sin(A) == 0.0" << endl;
+// if (sin(B) == 0.0) cout << "sin(B) == 0.0" << endl;
+// if (sin(C) == 0.0) cout << "sin(C) == 0.0" << endl;
+
+ double aire_triangle_spherique = rayon * rayon * (a + b + c - M_PI);
- double aire_triangle_spherique = rayon * rayon * (a + b + c - pi);
+ double aire_sphere = 4.0 * M_PI * rayon * rayon;
- double aire_sphere = 4 * pi * rayon * rayon;
+ return ( (4.0 * 0.3333333 * M_PI * rayon * rayon * rayon) * (aire_triangle_spherique / aire_sphere) ) ;
- // attention division par zero !!!!
- return ( (4 * 0.3333333 * pi * rayon * rayon * rayon) * (aire_triangle_spherique / aire_sphere) ) ;
-
}
Modified: trunk/extra/SpherePadder/SpherePadder.hpp
===================================================================
--- trunk/extra/SpherePadder/SpherePadder.hpp 2009-02-26 08:35:54 UTC (rev 1695)
+++ trunk/extra/SpherePadder/SpherePadder.hpp 2009-02-26 12:47:38 UTC (rev 1696)
@@ -123,6 +123,7 @@
// High level pading functions
void pad_5 ();
void tetra_pad();
+ void densify();
// void insert_sphere(double x, double y, double z, double R);
// void densify ();
Modified: trunk/extra/SpherePadder/TetraMesh.cpp
===================================================================
--- trunk/extra/SpherePadder/TetraMesh.cpp 2009-02-26 08:35:54 UTC (rev 1695)
+++ trunk/extra/SpherePadder/TetraMesh.cpp 2009-02-26 12:47:38 UTC (rev 1696)
@@ -1,322 +1,325 @@
/*************************************************************************
-* Copyright (C) 2009 by Jean Francois Jerier *
-* jerier@xxxxxxxxxxxxxxx *
-* Copyright (C) 2009 by Vincent Richefeu *
-* vincent.richefeu@xxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
+ * Copyright (C) 2009 by Jean Francois Jerier *
+ * jerier@xxxxxxxxxxxxxxx *
+ * Copyright (C) 2009 by Vincent Richefeu *
+ * vincent.richefeu@xxxxxxxxxxx *
+ * *
+ * This program is free software; it is licensed under the terms of the *
+ * GNU General Public License v2 or later. See file LICENSE for details. *
+ *************************************************************************/
#include "TetraMesh.hpp"
TetraMesh::TetraMesh ()
{
- isOrganized = false;
+ isOrganized = false;
}
+
void TetraMesh::read_gmsh (const char* name)
{
- ifstream meshFile(name);
- if(!meshFile)
- {
- cerr << "TetraMesh::read_gmsh, cannot open file " << name << endl;
- return;
- }
-
- string token;
- char not_read[150];
- meshFile >> token;
-
- while(meshFile)
- {
- if (token == "$Nodes")
+ ifstream meshFile(name);
+ if(!meshFile)
{
- unsigned int nbnodes;
- unsigned int num_node;
- Node N;
-
- meshFile >> nbnodes;
- for (unsigned int n = 0 ; n < nbnodes ; ++n)
- {
- meshFile >> num_node >> N.x >> N.y >> N.z;
-cout << num_node << " " << N.x << "" << N.y << " " << N.z << endl;
- node.push_back(N);
- }
+ cerr << "TetraMesh::read_gmsh, cannot open file " << name << endl;
+ return;
}
- if (token == "$Elements")
- {
- unsigned int nbElements;
- unsigned int num_element, element_type, nbTags ;
- Tetraedre T;
- unsigned int t = 0;
+ string token;
+ char not_read[150];
+ meshFile >> token;
- meshFile >> nbElements;
- for (unsigned int e = 0 ; e < nbElements ; ++e)
- {
- meshFile >> num_element >> element_type;
- if (element_type != 4) // 4-node tetrahedron
+ while(meshFile)
+ {
+ if (token == "$Nodes")
{
- meshFile.getline(not_read,150);
- continue;
+ unsigned int nbnodes;
+ unsigned int num_node;
+ Node N;
+
+ meshFile >> nbnodes;
+ for (unsigned int n = 0 ; n < nbnodes ; ++n)
+ {
+ meshFile >> num_node >> N.x >> N.y >> N.z;
+ // cout << num_node << " " << N.x << "" << N.y << " " << N.z << endl;
+ node.push_back(N);
+ }
}
+ if (token == "$Elements")
+ {
+ unsigned int nbElements;
+ unsigned int num_element, element_type, nbTags ;
+ Tetraedre T;
+ unsigned int t = 0;
- meshFile >> nbTags;
- // the third tag is the number of a mesh partition to which the element belongs
- unsigned int tag;
- for (unsigned int tg = 0 ; tg < nbTags ; ++(tg))
- { meshFile >> tag; }
+ meshFile >> nbElements;
+ for (unsigned int e = 0 ; e < nbElements ; ++e)
+ {
+ meshFile >> num_element >> element_type;
+ // 4-node tetrahedron
+ if (element_type != 4)
+ {
+ meshFile.getline(not_read,150);
+ continue;
+ }
- meshFile >> T.nodeId[0] >> T.nodeId[1] >> T.nodeId[2] >> T.nodeId[3];
+ meshFile >> nbTags;
+ // the third tag is the number of a mesh partition to which the element belongs
+ unsigned int tag;
+ for (unsigned int tg = 0 ; tg < nbTags ; ++(tg))
+ { meshFile >> tag; }
- // numbers begin at 0 instead of 1
- // (0 in C/C++ corresponds to 1 in the file)
- T.nodeId[0] -= 1;
- T.nodeId[1] -= 1;
- T.nodeId[2] -= 1;
- T.nodeId[3] -= 1;
+ meshFile >> T.nodeId[0] >> T.nodeId[1] >> T.nodeId[2] >> T.nodeId[3];
- node[T.nodeId[0]].tetraOwner.push_back(t);
- node[T.nodeId[1]].tetraOwner.push_back(t);
- node[T.nodeId[2]].tetraOwner.push_back(t);
- node[T.nodeId[3]].tetraOwner.push_back(t);
+ // numbers begin at 0 instead of 1
+ // (0 in C/C++ corresponds to 1 in the file)
+ T.nodeId[0] -= 1;
+ T.nodeId[1] -= 1;
+ T.nodeId[2] -= 1;
+ T.nodeId[3] -= 1;
- tetraedre.push_back(T);
- ++t;
- }
- }
+ node[T.nodeId[0]].tetraOwner.push_back(t);
+ node[T.nodeId[1]].tetraOwner.push_back(t);
+ node[T.nodeId[2]].tetraOwner.push_back(t);
+ node[T.nodeId[3]].tetraOwner.push_back(t);
- if (token == "$EndElements") break;
+ tetraedre.push_back(T);
+ ++t;
+ }
+ }
- meshFile >> token;
- }
+ if (token == "$EndElements") break;
- organize ();
+ meshFile >> token;
+ }
+
+ organize ();
}
+
void TetraMesh::read (const char* name)
{
- ifstream meshFile(name);
- if(!meshFile)
- {
- cerr << "TetraMesh::read_data, cannot open file " << name << endl;
- return;
- }
+ ifstream meshFile(name);
+ if(!meshFile)
+ {
+ cerr << "TetraMesh::read_data, cannot open file " << name << endl;
+ return;
+ }
- string token;
- meshFile >> token;
- bool with_numbers = false;
- unsigned int number;
+ string token;
+ meshFile >> token;
+ bool with_numbers = false;
+ unsigned int number;
- while(meshFile)
- {
- if (token == "WITH_NUMBERS") with_numbers = true;
- if (token == "NODES")
- {
- unsigned int nbnodes;
- Node N;
+ while(meshFile)
+ {
+ if (token == "WITH_NUMBERS") with_numbers = true;
+ if (token == "NODES")
+ {
+ unsigned int nbnodes;
+ Node N;
- meshFile >> nbnodes;
- for (unsigned int n = 0 ; n < nbnodes ; ++n)
- {
- if (with_numbers) meshFile >> number;
- meshFile >> N.x >> N.y >> N.z;
+ meshFile >> nbnodes;
+ for (unsigned int n = 0 ; n < nbnodes ; ++n)
+ {
+ if (with_numbers) meshFile >> number;
+ meshFile >> N.x >> N.y >> N.z;
- //cout << n << " " << N.x << " " << N.y << " " << N.z << endl;
- node.push_back(N);
- }
- }
+ //cout << n << " " << N.x << " " << N.y << " " << N.z << endl;
+ node.push_back(N);
+ }
+ }
- if (token == "TETRA")
- {
- unsigned int nbTetra;
-unsigned int nbignored = 0;
- Tetraedre T;
+ if (token == "TETRA")
+ {
+ unsigned int nbTetra;
+ unsigned int nbignored = 0;
+ Tetraedre T;
- meshFile >> nbTetra;
- for (unsigned int t = 0 ; t < nbTetra ; ++t)
- {
- if (with_numbers) meshFile >> number;
- meshFile >> T.nodeId[0] >> T.nodeId[1] >> T.nodeId[2] >> T.nodeId[3];
+ meshFile >> nbTetra;
+ for (unsigned int t = 0 ; t < nbTetra ; ++t)
+ {
+ if (with_numbers) meshFile >> number;
+ meshFile >> T.nodeId[0] >> T.nodeId[1] >> T.nodeId[2] >> T.nodeId[3];
- //cout << t << " | " << T.nodeId[0] << " " << T.nodeId[1] << " " << T.nodeId[2] << " " << T.nodeId[3] << endl;
+ //cout << t << " | " << T.nodeId[0] << " " << T.nodeId[1] << " " << T.nodeId[2] << " " << T.nodeId[3] << endl;
- // numbers begin at 0 instead of 1
- // (0 in C/C++ corresponds to 1 in the file)
- T.nodeId[0] -= 1;
- T.nodeId[1] -= 1;
- T.nodeId[2] -= 1;
- T.nodeId[3] -= 1;
+ // numbers begin at 0 instead of 1
+ // (0 in C/C++ corresponds to 1 in the file)
+ T.nodeId[0] -= 1;
+ T.nodeId[1] -= 1;
+ T.nodeId[2] -= 1;
+ T.nodeId[3] -= 1;
- if (T.nodeId[0] >= node.size() || T.nodeId[1] >= node.size() || T.nodeId[2] >= node.size() || T.nodeId[3] >= node.size())
- {
- nbignored++;
- //cerr << "(1) error in mesh file, tetrahedron ignored! nbignored = " << nbignored << endl;
- continue;
- }
- if (T.nodeId[0] < 0 || T.nodeId[1] < 0 || T.nodeId[2] < 0 || T.nodeId[3] < 0)
- {
- nbignored++;
- //cerr << "(2) error in mesh file, tetrahedron ignored! nbignored = " << nbignored << endl;
- continue;
- }
- //cout << "num = " << (t-nbignored) << endl;
+ if (T.nodeId[0] >= node.size() || T.nodeId[1] >= node.size() || T.nodeId[2] >= node.size() || T.nodeId[3] >= node.size())
+ {
+ nbignored++;
+ //cerr << "(1) error in mesh file, tetrahedron ignored! nbignored = " << nbignored << endl;
+ continue;
+ }
+ if (T.nodeId[0] < 0 || T.nodeId[1] < 0 || T.nodeId[2] < 0 || T.nodeId[3] < 0)
+ {
+ nbignored++;
+ //cerr << "(2) error in mesh file, tetrahedron ignored! nbignored = " << nbignored << endl;
+ continue;
+ }
+ //cout << "num = " << (t-nbignored) << endl;
- node[T.nodeId[0]].tetraOwner.push_back(t-nbignored);
- node[T.nodeId[1]].tetraOwner.push_back(t-nbignored);
- node[T.nodeId[2]].tetraOwner.push_back(t-nbignored);
- node[T.nodeId[3]].tetraOwner.push_back(t-nbignored);
+ node[T.nodeId[0]].tetraOwner.push_back(t-nbignored);
+ node[T.nodeId[1]].tetraOwner.push_back(t-nbignored);
+ node[T.nodeId[2]].tetraOwner.push_back(t-nbignored);
+ node[T.nodeId[3]].tetraOwner.push_back(t-nbignored);
- tetraedre.push_back(T);
+ tetraedre.push_back(T);
- }
- }
+ }
+ }
- if (token == "EOF") break;
+ if (token == "EOF") break;
- meshFile >> token;
- }
+ meshFile >> token;
+ }
- organize ();
+ organize ();
}
+
void TetraMesh::write_surface_MGP (const char* name)
{
- ofstream fmgpost(name);
+ ofstream fmgpost(name);
- fmgpost << "<?xml version=\"1.0\"?>" << endl
- << " <mgpost mode=\"3D\">" << endl
- << " <state id=\"" << 1
- << "\" time=\"" << 0.0 << "\">" << endl;
+ fmgpost << "<?xml version=\"1.0\"?>" << endl
+ << " <mgpost mode=\"3D\">" << endl
+ << " <state id=\"" << 1
+ << "\" time=\"" << 0.0 << "\">" << endl;
- double xm,ym,zm;
- unsigned int id1,id2,id3;
- unsigned int idPolyhedron = 1;
- double div3 = 0.3333333333333333;
- for (unsigned int i = 0 ; i < face.size() ; ++i)
- {
- if(!(face[i].belongBoundary)) continue;
- id1 = face[i].nodeId[0];
- id2 = face[i].nodeId[1];
- id3 = face[i].nodeId[2];
- xm = div3 * (node[id1].x + node[id2].x + node[id3].x);
- ym = div3 * (node[id1].y + node[id2].y + node[id3].y);
- zm = div3 * (node[id1].z + node[id2].z + node[id3].z);
+ double xm,ym,zm;
+ unsigned int id1,id2,id3;
+ unsigned int idPolyhedron = 1;
+ double div3 = 0.3333333333333333;
+ for (unsigned int i = 0 ; i < face.size() ; ++i)
+ {
+ if(!(face[i].belongBoundary)) continue;
+ id1 = face[i].nodeId[0];
+ id2 = face[i].nodeId[1];
+ id3 = face[i].nodeId[2];
+ xm = div3 * (node[id1].x + node[id2].x + node[id3].x);
+ ym = div3 * (node[id1].y + node[id2].y + node[id3].y);
+ zm = div3 * (node[id1].z + node[id2].z + node[id3].z);
- fmgpost << " <body>" << endl;
- fmgpost << " <POLYE id=\"" << idPolyhedron++ << "\" r=\"" << min_segment_length << "\">" << endl
+ fmgpost << " <body>" << endl;
+ fmgpost << " <POLYE id=\"" << idPolyhedron++ << "\" r=\"" << min_segment_length << "\">" << endl
<< " <position x=\"" << xm + xtrans << "\" y=\"" << ym + ytrans << "\" z=\"" << zm + ztrans << "\"/>" << endl
<< " <node x=\"" << node[id1].x - xm << "\" y=\"" << node[id1].y - ym << "\" z=\"" << node[id1].z - zm << "\"/>" << endl
<< " <node x=\"" << node[id2].x - xm << "\" y=\"" << node[id2].y - ym << "\" z=\"" << node[id2].z - zm << "\"/>" << endl
<< " <node x=\"" << node[id3].x - xm << "\" y=\"" << node[id3].y - ym << "\" z=\"" << node[id3].z - zm << "\"/>" << endl;
- if (face[i].normal_swap)
- fmgpost << " <face n1=\"" << 0 << "\" n2=\"" << 2 << "\" n3=\"" << 1 << "\"/>" << endl;
- else
- fmgpost << " <face n1=\"" << 0 << "\" n2=\"" << 1 << "\" n3=\"" << 2 << "\"/>" << endl;
+ if (face[i].normal_swap)
+ fmgpost << " <face n1=\"" << 0 << "\" n2=\"" << 2 << "\" n3=\"" << 1 << "\"/>" << endl;
+ else
+ fmgpost << " <face n1=\"" << 0 << "\" n2=\"" << 1 << "\" n3=\"" << 2 << "\"/>" << endl;
- fmgpost << " </POLYE>" << endl << flush;
+ fmgpost << " </POLYE>" << endl << flush;
- fmgpost << " </body>" << endl;
- }
-
- fmgpost << " </state>" << endl
- << " </mgpost>" << endl;
+ fmgpost << " </body>" << endl;
+ }
+
+ fmgpost << " </state>" << endl
+ << " </mgpost>" << endl;
}
-
int compareInt (const void * a, const void * b)
{
- return ( *(int*)a > *(int*)b ) ? 1 :-1;
+ return ( *(int*)a > *(int*)b ) ? 1 :-1;
}
+
void TetraMesh::organize ()
{
- // Translate all nodes in such a manner that all coordinates are > 0
- xtrans = node[0].x;
- ytrans = node[0].y;
- ztrans = node[0].z;
- for (unsigned int i = 1 ; i < node.size() ; ++i)
- {
- xtrans = (xtrans < node[i].x) ? xtrans : node[i].x;
- ytrans = (ytrans < node[i].y) ? ytrans : node[i].y;
- ztrans = (ztrans < node[i].z) ? ztrans : node[i].z;
- }
- xtrans = (xtrans < 0.0) ? -xtrans : 0.0;
- ytrans = (ytrans < 0.0) ? -ytrans : 0.0;
- ztrans = (ztrans < 0.0) ? -ztrans : 0.0;
- for (unsigned int i = 0 ; i < node.size() ; ++i)
- {
- node[i].x += xtrans;
- node[i].y += ytrans;
- node[i].z += ztrans;
- }
+ // Translate all nodes in such a manner that all coordinates are > 0
+ xtrans = node[0].x;
+ ytrans = node[0].y;
+ ztrans = node[0].z;
+ for (unsigned int i = 1 ; i < node.size() ; ++i)
+ {
+ xtrans = (xtrans < node[i].x) ? xtrans : node[i].x;
+ ytrans = (ytrans < node[i].y) ? ytrans : node[i].y;
+ ztrans = (ztrans < node[i].z) ? ztrans : node[i].z;
+ }
+ xtrans = (xtrans < 0.0) ? -xtrans : 0.0;
+ ytrans = (ytrans < 0.0) ? -ytrans : 0.0;
+ ztrans = (ztrans < 0.0) ? -ztrans : 0.0;
+ for (unsigned int i = 0 ; i < node.size() ; ++i)
+ {
+ node[i].x += xtrans;
+ node[i].y += ytrans;
+ node[i].z += ztrans;
+ }
- // Organize tetraedre nodes in ascending order
- for (unsigned int i = 0 ; i < tetraedre.size() ; ++i)
- {
- qsort (&(tetraedre[i].nodeId[0]), 4, sizeof(int), compareInt);
- }
+ // Organize tetraedre nodes in ascending order
+ for (unsigned int i = 0 ; i < tetraedre.size() ; ++i)
+ {
+ qsort (&(tetraedre[i].nodeId[0]), 4, sizeof(int), compareInt);
+ }
- // Face creation
- vector <Face> tmpFace; // This will contain all faces more than one time (with duplications)
- Face F;
- F.tetraOwner.push_back(0);
- F.belongBoundary = true;
+ // Face creation
+ vector <Face> tmpFace; // This will contain all faces more than one time (with duplications)
+ Face F;
+ F.tetraOwner.push_back(0);
+ F.belongBoundary = true;
F.normal_swap = false;
- for (unsigned int i = 0 ; i < tetraedre.size() ; ++i)
- {
- F.tetraOwner[0] = i;
+ for (unsigned int i = 0 ; i < tetraedre.size() ; ++i)
+ {
+ F.tetraOwner[0] = i;
- // Notice that nodes are still organized in ascending order
- F.nodeId[0] = tetraedre[i].nodeId[0];
- F.nodeId[1] = tetraedre[i].nodeId[1];
- F.nodeId[2] = tetraedre[i].nodeId[2];
- tmpFace.push_back(F);
+ // Notice that nodes are still organized in ascending order
+ F.nodeId[0] = tetraedre[i].nodeId[0];
+ F.nodeId[1] = tetraedre[i].nodeId[1];
+ F.nodeId[2] = tetraedre[i].nodeId[2];
+ tmpFace.push_back(F);
- F.nodeId[0] = tetraedre[i].nodeId[1];
- F.nodeId[1] = tetraedre[i].nodeId[2];
- F.nodeId[2] = tetraedre[i].nodeId[3];
- tmpFace.push_back(F);
+ F.nodeId[0] = tetraedre[i].nodeId[1];
+ F.nodeId[1] = tetraedre[i].nodeId[2];
+ F.nodeId[2] = tetraedre[i].nodeId[3];
+ tmpFace.push_back(F);
- F.nodeId[0] = tetraedre[i].nodeId[0];
- F.nodeId[1] = tetraedre[i].nodeId[2];
- F.nodeId[2] = tetraedre[i].nodeId[3];
- tmpFace.push_back(F);
+ F.nodeId[0] = tetraedre[i].nodeId[0];
+ F.nodeId[1] = tetraedre[i].nodeId[2];
+ F.nodeId[2] = tetraedre[i].nodeId[3];
+ tmpFace.push_back(F);
- F.nodeId[0] = tetraedre[i].nodeId[0];
- F.nodeId[1] = tetraedre[i].nodeId[1];
- F.nodeId[2] = tetraedre[i].nodeId[3];
- tmpFace.push_back(F);
- }
+ F.nodeId[0] = tetraedre[i].nodeId[0];
+ F.nodeId[1] = tetraedre[i].nodeId[1];
+ F.nodeId[2] = tetraedre[i].nodeId[3];
+ tmpFace.push_back(F);
+ }
- face.push_back(tmpFace[0]);
- bool duplicated;
- for (unsigned int i = 1 ; i < tmpFace.size() ; ++i)
- {
- duplicated = false;
- for (unsigned int n = 0 ; n < face.size() ; ++n)
- {
- if ( tmpFace[i].nodeId[0] == face[n].nodeId[0]
- && tmpFace[i].nodeId[1] == face[n].nodeId[1]
- && tmpFace[i].nodeId[2] == face[n].nodeId[2])
- {
- duplicated = true;
- face[n].tetraOwner.push_back(tmpFace[i].tetraOwner[0]);
- face[n].belongBoundary = false;
- break;
- }
- }
+ face.push_back(tmpFace[0]);
+ bool duplicated;
+ for (unsigned int i = 1 ; i < tmpFace.size() ; ++i)
+ {
+ duplicated = false;
+ for (unsigned int n = 0 ; n < face.size() ; ++n)
+ {
+ if ( tmpFace[i].nodeId[0] == face[n].nodeId[0]
+ && tmpFace[i].nodeId[1] == face[n].nodeId[1]
+ && tmpFace[i].nodeId[2] == face[n].nodeId[2])
+ {
+ duplicated = true;
+ face[n].tetraOwner.push_back(tmpFace[i].tetraOwner[0]);
+ face[n].belongBoundary = false;
+ break;
+ }
+ }
- if (!duplicated)
- {
- face.push_back(tmpFace[i]);
- }
- }
- tmpFace.clear(); // It should be usefull to free memory before the end of the function
+ if (!duplicated)
+ {
+ face.push_back(tmpFace[i]);
+ }
+ }
+ tmpFace.clear(); // It should be usefull to free memory before the end of the function
// Definition of unit vectors normal to the boundary faces
unsigned int n1,n2,n3,n4;
@@ -326,186 +329,181 @@
double scal_product;
for (unsigned int f = 0 ; f < face.size() ; ++f)
{
- if (!(face[f].belongBoundary)) continue;
+ if (!(face[f].belongBoundary)) continue;
- current_tetra = tetraedre[ face[f].tetraOwner[0] ];
+ current_tetra = tetraedre[ face[f].tetraOwner[0] ];
- n1 = current_tetra.nodeId[0];
- n2 = current_tetra.nodeId[1];
- n3 = current_tetra.nodeId[2];
- n4 = current_tetra.nodeId[3];
+ n1 = current_tetra.nodeId[0];
+ n2 = current_tetra.nodeId[1];
+ n3 = current_tetra.nodeId[2];
+ n4 = current_tetra.nodeId[3];
- s1 = face[f].nodeId[0];
- s2 = face[f].nodeId[1];
- s3 = face[f].nodeId[2];
+ s1 = face[f].nodeId[0];
+ s2 = face[f].nodeId[1];
+ s3 = face[f].nodeId[2];
- if ((n1 != s1) && (n1 != s2) && (n1 != s3)) {s4 = n1;}
- else if ((n2 != s1) && (n2 != s2) && (n2 != s3)) {s4 = n2;}
- else if ((n3 != s1) && (n3 != s2) && (n3 != s3)) {s4 = n3;}
- else {s4 = n4;}
+ if ((n1 != s1) && (n1 != s2) && (n1 != s3)) {s4 = n1;}
+ else if ((n2 != s1) && (n2 != s2) && (n2 != s3)) {s4 = n2;}
+ else if ((n3 != s1) && (n3 != s2) && (n3 != s3)) {s4 = n3;}
+ else {s4 = n4;}
- v12[0] = node[s2].x - node[s1].x;
- v12[1] = node[s2].y - node[s1].y;
- v12[2] = node[s2].z - node[s1].z;
+ v12[0] = node[s2].x - node[s1].x;
+ v12[1] = node[s2].y - node[s1].y;
+ v12[2] = node[s2].z - node[s1].z;
- v13[0] = node[s3].x - node[s1].x;
- v13[1] = node[s3].y - node[s1].y;
- v13[2] = node[s3].z - node[s1].z;
+ v13[0] = node[s3].x - node[s1].x;
+ v13[1] = node[s3].y - node[s1].y;
+ v13[2] = node[s3].z - node[s1].z;
- v14[0] = node[s4].x - node[s1].x;
- v14[1] = node[s4].y - node[s1].y;
- v14[2] = node[s4].z - node[s1].z;
+ v14[0] = node[s4].x - node[s1].x;
+ v14[1] = node[s4].y - node[s1].y;
+ v14[2] = node[s4].z - node[s1].z;
- vect_product[0] = v12[1]*v13[2] - v12[2]*v13[1];
- vect_product[1] = v12[2]*v13[0] - v12[0]*v13[2];
- vect_product[2] = v12[0]*v13[1] - v12[1]*v13[0];
+ vect_product[0] = v12[1]*v13[2] - v12[2]*v13[1];
+ vect_product[1] = v12[2]*v13[0] - v12[0]*v13[2];
+ vect_product[2] = v12[0]*v13[1] - v12[1]*v13[0];
- scal_product = v14[0]*vect_product[0] + v14[1]*vect_product[1] + v14[2]*vect_product[2];
+ scal_product = v14[0]*vect_product[0] + v14[1]*vect_product[1] + v14[2]*vect_product[2];
- if (scal_product > 0.0) face[f].normal_swap = true;
+ if (scal_product > 0.0) face[f].normal_swap = true;
}
-
// Define for each node the faces it belong to
- for (unsigned int f = 0 ; f < face.size() ; ++f)
- {
- node[face[f].nodeId[0]].faceOwner.push_back(f);
- node[face[f].nodeId[1]].faceOwner.push_back(f);
- node[face[f].nodeId[2]].faceOwner.push_back(f);
- }
+ for (unsigned int f = 0 ; f < face.size() ; ++f)
+ {
+ node[face[f].nodeId[0]].faceOwner.push_back(f);
+ node[face[f].nodeId[1]].faceOwner.push_back(f);
+ node[face[f].nodeId[2]].faceOwner.push_back(f);
+ }
- // Segment creation
- vector <Segment> tmpSegment;
- Segment S;
- S.faceOwner.push_back(0);
+ // Segment creation
+ vector <Segment> tmpSegment;
+ Segment S;
+ S.faceOwner.push_back(0);
- for (unsigned int i = 0 ; i < face.size() ; ++i)
- {
- S.faceOwner[0] = i;
- S.nodeId[0] = face[i].nodeId[0];
- S.nodeId[1] = face[i].nodeId[1];
- tmpSegment.push_back(S);
+ for (unsigned int i = 0 ; i < face.size() ; ++i)
+ {
+ S.faceOwner[0] = i;
+ S.nodeId[0] = face[i].nodeId[0];
+ S.nodeId[1] = face[i].nodeId[1];
+ tmpSegment.push_back(S);
- S.nodeId[0] = face[i].nodeId[1];
- S.nodeId[1] = face[i].nodeId[2];
- tmpSegment.push_back(S);
+ S.nodeId[0] = face[i].nodeId[1];
+ S.nodeId[1] = face[i].nodeId[2];
+ tmpSegment.push_back(S);
- S.nodeId[0] = face[i].nodeId[0];
- S.nodeId[1] = face[i].nodeId[2];
- tmpSegment.push_back(S);
- }
+ S.nodeId[0] = face[i].nodeId[0];
+ S.nodeId[1] = face[i].nodeId[2];
+ tmpSegment.push_back(S);
+ }
- segment.push_back(tmpSegment[0]);
- for (unsigned int i = 1 ; i < tmpSegment.size() ; ++i)
- {
- duplicated = false;
- for (unsigned int n = 0 ; n < segment.size() ; ++n)
- {
- if ( tmpSegment[i].nodeId[0] == segment[n].nodeId[0]
- && tmpSegment[i].nodeId[1] == segment[n].nodeId[1])
- {
- duplicated = true;
- segment[n].faceOwner.push_back(tmpSegment[i].faceOwner[0]);
- break;
- }
+ segment.push_back(tmpSegment[0]);
+ for (unsigned int i = 1 ; i < tmpSegment.size() ; ++i)
+ {
+ duplicated = false;
+ for (unsigned int n = 0 ; n < segment.size() ; ++n)
+ {
+ if ( tmpSegment[i].nodeId[0] == segment[n].nodeId[0]
+ && tmpSegment[i].nodeId[1] == segment[n].nodeId[1])
+ {
+ duplicated = true;
+ segment[n].faceOwner.push_back(tmpSegment[i].faceOwner[0]);
+ break;
+ }
- }
+ }
- if (!duplicated)
- {
- segment.push_back(tmpSegment[i]);
- }
- }
+ if (!duplicated)
+ {
+ segment.push_back(tmpSegment[i]);
+ }
+ }
- for (unsigned int s = 0 ; s < segment.size() ; ++s)
- {
- node[segment[s].nodeId[0]].segmentOwner.push_back(s);
- node[segment[s].nodeId[1]].segmentOwner.push_back(s);
- }
+ for (unsigned int s = 0 ; s < segment.size() ; ++s)
+ {
+ node[segment[s].nodeId[0]].segmentOwner.push_back(s);
+ node[segment[s].nodeId[1]].segmentOwner.push_back(s);
+ }
- // Compute length of segments
- double lx,ly,lz;
- unsigned int id1,id2;
+ // Compute length of segments
+ double lx,ly,lz;
+ unsigned int id1,id2;
- // length_min = length(0)
- id1 = segment[0].nodeId[0];
- id2 = segment[0].nodeId[1];
+ // length_min = length(0)
+ id1 = segment[0].nodeId[0];
+ id2 = segment[0].nodeId[1];
+ lx = node[id2].x - node[id1].x;
+ ly = node[id2].y - node[id1].y;
+ lz = node[id2].z - node[id1].z;
+ min_segment_length = sqrt(lx*lx + ly*ly + lz*lz);
+
+ mean_segment_length = 0.0;
+ max_segment_length = 0.0;
+
+ for (unsigned int s = 0 ; s < segment.size() ; ++s)
+ {
+ id1 = segment[s].nodeId[0];
+ id2 = segment[s].nodeId[1];
+
lx = node[id2].x - node[id1].x;
ly = node[id2].y - node[id1].y;
lz = node[id2].z - node[id1].z;
- min_segment_length = sqrt(lx*lx + ly*ly + lz*lz);
+ segment[s].length = sqrt(lx*lx + ly*ly + lz*lz);
- mean_segment_length = 0.0;
- max_segment_length = 0.0;
+ mean_segment_length += segment[s].length;
+ min_segment_length = (segment[s].length < min_segment_length) ? segment[s].length : min_segment_length;
+ max_segment_length = (segment[s].length > max_segment_length) ? segment[s].length : max_segment_length;
+ }
+ mean_segment_length /= (double)(segment.size());
- for (unsigned int s = 0 ; s < segment.size() ; ++s)
- {
- id1 = segment[s].nodeId[0];
- id2 = segment[s].nodeId[1];
+ cerr << "mean_segment_length = " << mean_segment_length << endl;
+ cerr << "min_segment_length = " << min_segment_length << endl;
+ cerr << "max_segment_length = " << max_segment_length << endl;
- lx = node[id2].x - node[id1].x;
- ly = node[id2].y - node[id1].y;
- lz = node[id2].z - node[id1].z;
-
- segment[s].length = sqrt(lx*lx + ly*ly + lz*lz);
-
- mean_segment_length += segment[s].length;
- min_segment_length = (segment[s].length < min_segment_length) ? segment[s].length : min_segment_length;
- max_segment_length = (segment[s].length > max_segment_length) ? segment[s].length : max_segment_length;
- }
- mean_segment_length /= (double)(segment.size());
-
- cerr << "mean_segment_length = " << mean_segment_length << endl;
- cerr << "min_segment_length = " << min_segment_length << endl;
- cerr << "max_segment_length = " << max_segment_length << endl;
-
- // Define tetraedre neighbors
- bool stop = false;
- for (unsigned int t1 = 0 ; t1 < tetraedre.size() ; ++t1)
- {
- tetraedre[t1].tetraNeighbor.push_back(t1);
- for (unsigned int t2 = 0 ; t2 < tetraedre.size() ; ++t2)
- {
- /*
- if ( (tetraedre[t1].nodeId[0] > tetraedre[t2].nodeId[3])
- || (tetraedre[t1].nodeId[3] < tetraedre[t2].nodeId[0]) ) continue;
- */
-
- // FIXME mettre du while... (?)
- stop = false;
- for (unsigned int i = 0 ; i < 4 ; i++)
- {
- for (unsigned int j = 0 ; j < 4 ; j++)
- {
- if (tetraedre[t1].nodeId[i] == tetraedre[t2].nodeId[j])
- {
- if (t1 != t2) tetraedre[t1].tetraNeighbor.push_back(t2);
- //if (t1 != t2) tetraedre[t2].tetraNeighbor.push_back(t1);
- stop = true;
- break;
- }
- }
- if (stop) break;
- }
- }
- }
-
-/*
- for (unsigned int t1 = 0 ; t1 < tetraedre.size() ; ++t1)
+ // Define tetraedre neighbors
+ bool stop = false;
+ for (unsigned int t1 = 0 ; t1 < tetraedre.size() ; ++t1)
+ {
+ tetraedre[t1].tetraNeighbor.push_back(t1);
+ for (unsigned int t2 = 0 ; t2 < tetraedre.size() ; ++t2)
{
- cerr << "Tetra " << t1 << " a " << tetraedre[t1].tetraNeighbor.size() << " voisins" << endl;
- for (unsigned int i = 0 ; i < tetraedre[t1].tetraNeighbor.size() ; i++)
- {
- cerr << tetraedre[t1].tetraNeighbor[i] << ' ';
- }
- cerr << endl;
+ /*
+ if ( (tetraedre[t1].nodeId[0] > tetraedre[t2].nodeId[3])
+ || (tetraedre[t1].nodeId[3] < tetraedre[t2].nodeId[0]) ) continue;
+ */
+
+ // FIXME mettre du while... (?)
+ stop = false;
+ for (unsigned int i = 0 ; i < 4 ; i++)
+ {
+ for (unsigned int j = 0 ; j < 4 ; j++)
+ {
+ if (tetraedre[t1].nodeId[i] == tetraedre[t2].nodeId[j])
+ {
+ if (t1 != t2) tetraedre[t1].tetraNeighbor.push_back(t2);
+ //if (t1 != t2) tetraedre[t2].tetraNeighbor.push_back(t1);
+ stop = true;
+ break;
+ }
+ }
+ if (stop) break;
+ }
}
- exit(0);
-*/
+ }
- isOrganized = true;
-}
+ /*
+ for (unsigned int t1 = 0 ; t1 < tetraedre.size() ; ++t1)
+ {
+ cerr << "Tetra " << t1 << " a " << tetraedre[t1].tetraNeighbor.size() << " voisins" << endl;
+ for (unsigned int i = 0 ; i < tetraedre[t1].tetraNeighbor.size() ; i++)
+ {
+ cerr << tetraedre[t1].tetraNeighbor[i] << ' ';
+ }
+ cerr << endl;
+ }
+ exit(0);
+ */
-
-
+ isOrganized = true;
+}
Modified: trunk/extra/SpherePadder/main.cpp
===================================================================
--- trunk/extra/SpherePadder/main.cpp 2009-02-26 08:35:54 UTC (rev 1695)
+++ trunk/extra/SpherePadder/main.cpp 2009-02-26 12:47:38 UTC (rev 1696)
@@ -24,16 +24,16 @@
TetraMesh * mesh = new TetraMesh();
- //mesh->read_gmsh("meshes/test2.msh");
- mesh->read("meshes/tomo.tetra");
- //mesh->write_surface_MGP ("tomo.mgp");
+ //mesh->read_gmsh("meshes/cube1194.msh");
+ mesh->read("meshes/test.tetra");
+ //mesh->write_surface_MGP ("cube.mgp");
SpherePadder * padder = new SpherePadder();
padder->plugTetraMesh(mesh);
//padder->add_spherical_probe(0.7);
padder->pad_5();
- //padder->tetra_pad();
+ padder->densify();
padder->save_mgpost("mgp.out.001");
// padder->save_Rxyz("spheres.Rxyz");
