yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #00862
[svn] r1619 - trunk/extra/SpherePadder
Author: richefeu
Date: 2009-01-09 15:44:43 +0100 (Fri, 09 Jan 2009)
New Revision: 1619
Modified:
trunk/extra/SpherePadder/SpherePadder.cpp
trunk/extra/SpherePadder/SpherePadder.hpp
trunk/extra/SpherePadder/TetraMesh.cpp
trunk/extra/SpherePadder/TetraMesh.hpp
trunk/extra/SpherePadder/main.cpp
Log:
- SpherePadder update
- add a basic 'user friend' interface for sphere packing generation. The module can of course be used without this interface.
Modified: trunk/extra/SpherePadder/SpherePadder.cpp
===================================================================
--- trunk/extra/SpherePadder/SpherePadder.cpp 2009-01-09 13:09:26 UTC (rev 1618)
+++ trunk/extra/SpherePadder/SpherePadder.cpp 2009-01-09 14:44:43 UTC (rev 1619)
@@ -35,12 +35,10 @@
combination.push_back(lst);
}
- max_overlap_rate = 1e-7;
+ max_overlap_rate = 1e-3;
n1 = n2 = n3 = n4 = n5 = n_densify = 0;
-
- //FIXME :
- // si on utilise la class dans un autre programme c'est peu etre pas util
-
+ trace_functions = true;
+ meshIsPlugged = false;
//rmin = 1e-2;
//rmax = 2e-2;
//rmoy = 0.5 * (rmin + rmax);
@@ -55,268 +53,290 @@
void SpherePadder::plugTetraMesh (TetraMesh * pluggedMesh)
{
- mesh = pluggedMesh;
+ mesh = pluggedMesh;
+ meshIsPlugged = true;
- // TODO mettre ce qui suite dans une fonction 'init()'
- // Si l'utilisateur n'a choisi qu'une valeur de ratio,
- // on cree les valeur de rmin et rmax:
- if (rmoy == 0 && ratio > 0)
- {
- rmoy = 0.125 * mesh->mean_segment_length; // 1/8
- rmin = (2.0 * rmoy) / (ratio + 1.0);
- //rmax = (2.0 * ratio * rmoy) / (ratio + 1.0);
- rmax = 2.0 * rmoy - rmin;
- //dr = 0.5 * (rmax - rmin);
- dr = rmax - rmoy;
- }
+ // TODO mettre ce qui suite dans une fonction 'init()'
+ // Si l'utilisateur n'a choisi qu'une valeur de ratio,
+ // on cree les valeur de rmin et rmax:
+ if (rmoy == 0 && ratio > 0)
+ {
+ rmoy = 0.125 * mesh->mean_segment_length; // 1/8
+ rmin = (2.0 * rmoy) / (ratio + 1.0);
+ //rmax = (2.0 * ratio * rmoy) / (ratio + 1.0);
+ rmax = 2.0 * rmoy - rmin;
+ //dr = 0.5 * (rmax - rmin);
+ dr = rmax - rmoy;
+ }
}
void SpherePadder::pad_5 ()
{
- // TODO check if all si ok (mesh exist...)
- place_at_nodes();
- place_at_segment_middle();
- cancel_overlap();
- place_at_faces();
- place_at_tetra_centers();
- //place_at_tetra_vertexes ();
+ if (mesh == 0)
+ {
+ cerr << "SpherePadder::pad_5, no mesh defined!" << endl;
+ return;
+ }
+
+ if (!(mesh->isOrganized))
+ {
+ cerr << "SpherePadder::pad_5, mesh is not valid!" << endl;
+ return;
+ }
+
+ place_at_nodes();
+ place_at_segment_middle();
+ cancel_overlaps();
+ place_at_faces();
+ place_at_tetra_centers();
+ place_at_tetra_vertexes ();
- cerr << "nb spheres = " << sphere.size() << endl;
-
+ cerr << "Total number of spheres = " << sphere.size() << endl;
+
}
void SpherePadder::save_mgpost (const char* name)
{
- cerr << "save mgp... ";
+ BEGIN_FUNCTION ("Save mgp");
ofstream fmgpost(name);
-
+
+ double xtrans = mesh->xtrans;
+ double ytrans = mesh->ytrans;
+ double ztrans = mesh->ztrans;
+
fmgpost << "<?xml version=\"1.0\"?>" << endl
- << " <mgpost mode=\"3D\">" << endl
- << " <newcolor name=\"at nodes\"/>" << endl
- << " <newcolor name=\"at segments\"/>" << endl
- << " <newcolor name=\"at faces\"/>" << endl
- << " <newcolor name=\"at tetra centers\"/>" << endl
- << " <newcolor name=\"at tetra vertexes\"/>" << endl
- << " <state id=\"" << 1
- << "\" time=\"" << 0.0 << "\">" << endl;
+ << " <mgpost mode=\"3D\">" << endl
+ << " <newcolor name=\"at nodes\"/>" << endl
+ << " <newcolor name=\"at segments\"/>" << endl
+ << " <newcolor name=\"at faces\"/>" << endl
+ << " <newcolor name=\"at tetra centers\"/>" << endl
+ << " <newcolor name=\"at tetra vertexes\"/>" << endl
+ << " <state id=\"" << 1
+ << "\" time=\"" << 0.0 << "\">" << endl;
for (unsigned int i = 0 ; i < sphere.size() ; ++i)
{
- fmgpost << " <body>" << endl;
- fmgpost << " <SPHER id=\"" << i+1 << "\" col=\"" << sphere[i].type << "\" r=\"" << sphere[i].R << "\">" << endl
- << " <position x=\"" << sphere[i].x << "\" y=\"" << sphere[i].y << "\" z=\"" << sphere[i].z << "\"/>" << endl
- << " </SPHER>" << endl << flush;
+ fmgpost << " <body>" << endl;
+ fmgpost << " <SPHER id=\"" << i+1 << "\" col=\"" << sphere[i].type << "\" r=\"" << sphere[i].R << "\">" << endl
+ << " <position x=\"" << sphere[i].x + xtrans << "\" y=\""
+ << sphere[i].y + ytrans << "\" z=\"" << sphere[i].z + ztrans << "\"/>" << endl
+ << " </SPHER>" << endl << flush;
// tmp (bricolage)
- if (i < mesh->node.size())
- for (unsigned int s = 0 ; s < mesh->segment.size() ; ++s)
- {
- if (mesh->segment[s].nodeId[0] == i)
- {
- fmgpost << " <SPSPx antac=\"" << mesh->segment[s].nodeId[1] + 1 << "\"/>" << endl;
- }
- }
+ if (i < mesh->node.size())
+ for (unsigned int s = 0 ; s < mesh->segment.size() ; ++s)
+ {
+ if (mesh->segment[s].nodeId[0] == i)
+ {
+ fmgpost << " <SPSPx antac=\"" << mesh->segment[s].nodeId[1] + 1 << "\"/>" << endl;
+ }
+ }
- fmgpost << " </body>" << endl;
- }
+ fmgpost << " </body>" << endl;
+ }
fmgpost << " </state>" << endl
- << " </mgpost>" << endl;
-
- cerr << "Done" << endl;
+ << " </mgpost>" << endl;
+
+ END_FUNCTION;
}
void SpherePadder::save_Rxyz (const char* name)
{
- cerr << "save Rxyz... ";
+ BEGIN_FUNCTION("Save Rxyz");
ofstream file(name);
+
+ double xtrans = mesh->xtrans;
+ double ytrans = mesh->ytrans;
+ double ztrans = mesh->ztrans;
for (unsigned int i = 0 ; i < sphere.size() ; ++i)
{
- file << sphere[i].R << " " << sphere[i].x << " " << sphere[i].y << " " << sphere[i].z << endl;
+ file << sphere[i].R << " " << sphere[i].x + xtrans << " " << sphere[i].y + ytrans << " " << sphere[i].z + ztrans << endl;
}
- cerr << "Done" << endl;
+ END_FUNCTION;
}
void SpherePadder::place_at_nodes ()
{
- unsigned int segId;
- Sphere S;
- S.type = AT_NODE;
-
- cerr << "place at nodes... ";
-
- for (unsigned int n = 0 ; n < mesh->node.size() ; ++n)
- {
- S.x = mesh->node[n].x;
- S.y = mesh->node[n].y;
- S.z = mesh->node[n].z;
- S.R = mesh->segment[ mesh->node[n].segmentOwner[0] ].length;
- for (unsigned int i = 1 ; i < mesh->node[n].segmentOwner.size() ; ++i)
- {
- segId = mesh->node[n].segmentOwner[i];
- S.R = (S.R < mesh->segment[segId].length) ? S.R : mesh->segment[segId].length;
- }
- S.R /= 4.0;
-
- S.tetraOwner = mesh->node[n].tetraOwner[0];
- mesh->tetraedre[S.tetraOwner].sphereId.push_back(n);
-
- sphere.push_back(S); ++(n1);
- }
- cerr << "Done" << endl;
+ BEGIN_FUNCTION("Place at nodes");
+
+ unsigned int segId;
+ Sphere S;
+ S.type = AT_NODE;
+
+ for (unsigned int n = 0 ; n < mesh->node.size() ; ++n)
+ {
+ S.x = mesh->node[n].x;
+ S.y = mesh->node[n].y;
+ S.z = mesh->node[n].z;
+ S.R = mesh->segment[ mesh->node[n].segmentOwner[0] ].length;
+ for (unsigned int i = 1 ; i < mesh->node[n].segmentOwner.size() ; ++i)
+ {
+ segId = mesh->node[n].segmentOwner[i];
+ S.R = (S.R < mesh->segment[segId].length) ? S.R : mesh->segment[segId].length;
+ }
+ S.R /= 4.0;
+
+ S.tetraOwner = mesh->node[n].tetraOwner[0];
+ mesh->tetraedre[S.tetraOwner].sphereId.push_back(n);
+
+ sphere.push_back(S); ++(n1);
+ }
+
+ END_FUNCTION;
}
void SpherePadder::place_at_segment_middle ()
-{
- cerr << "place at segment middle... ";
- Sphere S;
- S.type = AT_SEGMENT;
- double x1,y1,z1;
- double x2,y2,z2;
- unsigned int id1,id2;
- unsigned int n0 = sphere.size();
-
- for (unsigned int s = 0 ; s < mesh->segment.size() ; ++s)
- {
- id1 = mesh->segment[s].nodeId[0];
- id2 = mesh->segment[s].nodeId[1];
-
- x1 = mesh->node[id1].x;
- y1 = mesh->node[id1].y;
- z1 = mesh->node[id1].z;
-
- x2 = mesh->node[id2].x;
- y2 = mesh->node[id2].y;
- z2 = mesh->node[id2].z;
-
- S.x = 0.5 * (x1 + x2);
- S.y = 0.5 * (y1 + y2);
- S.z = 0.5 * (z1 + z2);
- S.R = 0.125 * mesh->segment[s].length;
- if (S.R < rmin) S.R = rmin;
- else if (S.R > rmax) S.R = rmoy + dr * (double)rand()/(double)RAND_MAX;
-
- S.tetraOwner = mesh->node[id1].tetraOwner[0];
- mesh->tetraedre[S.tetraOwner].sphereId.push_back(n0 + s);
- sphere.push_back(S); ++(n2);
+{
+ BEGIN_FUNCTION("Place at segment middle");
+ Sphere S;
+ S.type = AT_SEGMENT;
+ double x1,y1,z1;
+ double x2,y2,z2;
+ unsigned int id1,id2;
+ unsigned int n0 = sphere.size();
+
+ for (unsigned int s = 0 ; s < mesh->segment.size() ; ++s)
+ {
+ id1 = mesh->segment[s].nodeId[0];
+ id2 = mesh->segment[s].nodeId[1];
- mesh->segment[s].sphereId = n0 + s;
-
- }
- cerr << "Done" << endl;
+ x1 = mesh->node[id1].x;
+ y1 = mesh->node[id1].y;
+ z1 = mesh->node[id1].z;
+
+ x2 = mesh->node[id2].x;
+ y2 = mesh->node[id2].y;
+ z2 = mesh->node[id2].z;
+
+ S.x = 0.5 * (x1 + x2);
+ S.y = 0.5 * (y1 + y2);
+ S.z = 0.5 * (z1 + z2);
+ S.R = 0.125 * mesh->segment[s].length;
+ if (S.R < rmin) S.R = rmin;
+ else if (S.R > rmax) S.R = rmoy + dr * (double)rand()/(double)RAND_MAX;
+
+ S.tetraOwner = mesh->node[id1].tetraOwner[0];
+ mesh->tetraedre[S.tetraOwner].sphereId.push_back(n0 + s);
+ sphere.push_back(S); ++(n2);
+
+ mesh->segment[s].sphereId = n0 + s;
+
+ }
+ END_FUNCTION;
}
void SpherePadder::place_at_faces ()
{
- cerr << "place at barycentre 3... ";
+ BEGIN_FUNCTION("Place at faces");
- // FIXME move the following loops in TetraMesh or at the end of place_at_segment_middle ??
- unsigned int sphereId,faceId;
- for (unsigned int s = 0 ; s < mesh->segment.size() ; ++s)
- {
- sphereId = mesh->segment[s].sphereId;
- for (unsigned int f = 0 ; f < mesh->segment[s].faceOwner.size() ; ++f)
- {
- faceId = mesh->segment[s].faceOwner[f];
- mesh->face[ faceId ].sphereId.push_back( sphereId );
- }
- }
-
- Sphere S;
- S.type = AT_FACE;
+ // FIXME move the following loops in TetraMesh or at the end of place_at_segment_middle ??
+ unsigned int sphereId,faceId;
+ for (unsigned int s = 0 ; s < mesh->segment.size() ; ++s)
+ {
+ sphereId = mesh->segment[s].sphereId;
+ for (unsigned int f = 0 ; f < mesh->segment[s].faceOwner.size() ; ++f)
+ {
+ faceId = mesh->segment[s].faceOwner[f];
+ mesh->face[ faceId ].sphereId.push_back( sphereId );
+ }
+ }
+
+ Sphere S;
+ S.type = AT_FACE;
- unsigned int ns = sphere.size();
- const double div3 = 0.3333333333333;
- Sphere S1,S2,S3;
+ unsigned int ns = sphere.size();
+ const double div3 = 0.3333333333333;
+ Sphere S1,S2,S3;
- for (unsigned int f = 0 ; f < mesh->face.size() ; ++f)
- {
+ for (unsigned int f = 0 ; f < mesh->face.size() ; ++f)
+ {
- S1 = sphere[ mesh->face[f].sphereId[0] ];
- S2 = sphere[ mesh->face[f].sphereId[1] ];
- S3 = sphere[ mesh->face[f].sphereId[2] ];
+ S1 = sphere[ mesh->face[f].sphereId[0] ];
+ S2 = sphere[ mesh->face[f].sphereId[1] ];
+ S3 = sphere[ mesh->face[f].sphereId[2] ];
- S.x = div3 * (S1.x + S2.x + S3.x);
- S.y = div3 * (S1.y + S2.y + S3.y);
- S.z = div3 * (S1.z + S2.z + S3.z);
- S.R = rmin;
-
- S.tetraOwner = mesh->node[ mesh->face[f].nodeId[0] ].tetraOwner[0];
- mesh->tetraedre[S.tetraOwner].sphereId.push_back(ns++);
- sphere.push_back(S); ++(n3);
- }
-
- for (unsigned int n = (n1+n2) ; n < sphere.size() ; ++n)
- {
- place_sphere_4contacts(n);
- }
+ S.x = div3 * (S1.x + S2.x + S3.x);
+ S.y = div3 * (S1.y + S2.y + S3.y);
+ S.z = div3 * (S1.z + S2.z + S3.z);
+ S.R = rmin;
+
+ S.tetraOwner = mesh->node[ mesh->face[f].nodeId[0] ].tetraOwner[0];
+ mesh->tetraedre[S.tetraOwner].sphereId.push_back(ns++);
+ sphere.push_back(S); ++(n3);
+ }
- cerr << "Done" << endl;
+ for (unsigned int n = (n1+n2) ; n < sphere.size() ; ++n)
+ {
+ place_sphere_4contacts(n);
+ }
+
+ END_FUNCTION;
}
double SpherePadder::distance_spheres(unsigned int i, unsigned int j)
{
- double lx,ly,lz;
- lx = sphere[j].x - sphere[i].x;
- ly = sphere[j].y - sphere[i].y;
- lz = sphere[j].z - sphere[i].z;
- return (sqrt(lx*lx + ly*ly + lz*lz) - sphere[i].R - sphere[j].R);
+ double lx,ly,lz;
+ lx = sphere[j].x - sphere[i].x;
+ ly = sphere[j].y - sphere[i].y;
+ lz = sphere[j].z - sphere[i].z;
+ return (sqrt(lx*lx + ly*ly + lz*lz) - sphere[i].R - sphere[j].R);
}
double SpherePadder::distance_centre_spheres(Sphere& S1, Sphere& S2)
{
- double lx,ly,lz;
- lx = S2.x - S1.x;
- ly = S2.y - S1.y;
- lz = S2.z - S1.z;
- return (sqrt(lx*lx + ly*ly + lz*lz));
+ double lx,ly,lz;
+ lx = S2.x - S1.x;
+ ly = S2.y - S1.y;
+ lz = S2.z - S1.z;
+ return (sqrt(lx*lx + ly*ly + lz*lz));
}
-void SpherePadder::cancel_overlap() // FIXME rename cancel_overlaps
+void SpherePadder::cancel_overlaps()
{
-
- cerr << "cancel_overlaps... ";
- unsigned int current_tetra_id,tetra_neighbor_id,j;
- Tetraedre current_tetra, tetra_neighbor;
- double distance,k;
- double distance_max = -max_overlap_rate * rmax;
-
- for(unsigned int i = 0 ; i < sphere.size(); ++i)
- {
- if (sphere[i].R < 0.0) continue;
- current_tetra_id = sphere[i].tetraOwner;
- current_tetra = mesh->tetraedre[current_tetra_id];
-
- for (unsigned int t = 0 ; t < current_tetra.tetraNeighbor.size() ; ++t)
- {
- tetra_neighbor_id = current_tetra.tetraNeighbor[t];
- tetra_neighbor = mesh->tetraedre[tetra_neighbor_id];
- for (unsigned int n = 0 ; n < tetra_neighbor.sphereId.size() ; ++n)
- {
- j = tetra_neighbor.sphereId[n];
+
+ BEGIN_FUNCTION("Cancel_overlaps");
+ unsigned int current_tetra_id,tetra_neighbor_id,j;
+ Tetraedre current_tetra, tetra_neighbor;
+ double distance,k;
+ double distance_max = -max_overlap_rate * rmax;
+
+ for(unsigned int i = 0 ; i < sphere.size(); ++i)
+ {
+ if (sphere[i].R < 0.0) continue;
+ current_tetra_id = sphere[i].tetraOwner;
+ current_tetra = mesh->tetraedre[current_tetra_id];
+
+ for (unsigned int t = 0 ; t < current_tetra.tetraNeighbor.size() ; ++t)
+ {
+ tetra_neighbor_id = current_tetra.tetraNeighbor[t];
+ tetra_neighbor = mesh->tetraedre[tetra_neighbor_id];
+ for (unsigned int n = 0 ; n < tetra_neighbor.sphereId.size() ; ++n)
+ {
+ j = tetra_neighbor.sphereId[n];
- if (sphere[j].R < 0.0) continue;
- if (i < j)
- {
- while ( (distance = distance_spheres(i,j)) < distance_max )
- {
- k = 1.0 + distance / (sphere[i].R + sphere[j].R);
- sphere[i].R *= k;
- sphere[j].R *= k;
- }
- }
- }
- }
- }
- cerr << "Done" << endl;
+ if (sphere[j].R < 0.0) continue;
+ if (i < j)
+ {
+ while ( (distance = distance_spheres(i,j)) < distance_max )
+ {
+ k = 1.0 + distance / (sphere[i].R + sphere[j].R);
+ sphere[i].R *= k;
+ sphere[j].R *= k;
+ }
+ }
+ }
+ }
+ }
+ END_FUNCTION;
}
@@ -400,9 +420,11 @@
// if ((distance_centre_spheres(S,sphere[n]) - (S.R + sphere[n].R)) < -max_overlap_rate * rmin) { failure = 128; break; }
// }
+ double distance_max = -max_overlap_rate * rmin;
for (unsigned n = 0 ; n < neighbor.size() ; ++n)
{
- if ((distance_centre_spheres(S,sphere[neighbor[n].sphereId]) - (S.R + sphere[neighbor[n].sphereId].R)) < -max_overlap_rate * rmin) { failure = 128; break; }
+ if ((distance_centre_spheres(S,sphere[neighbor[n].sphereId]) - (S.R + sphere[neighbor[n].sphereId].R)) < distance_max)
+ { failure = 128; break; }
}
}
@@ -417,7 +439,7 @@
}
}
- // sphere[sphereId].R = 0.0; //debug
+
return 0;
}
@@ -537,8 +559,10 @@
}
else return fail_det;
+ // FIXME use not sqrt to speed up... (squared_distance_vector3)
+ // for the moment it is not critical
+
// Check interpenetration between spheres
-
double distance1 = distance_vector3 (centre,C1) - (R + R1);
double distance2 = distance_vector3 (centre,C2) - (R + R2);
double distance3 = distance_vector3 (centre,C3) - (R + R3);
@@ -552,7 +576,7 @@
{ return fail_overlap; }
// The gap between spheres must not be too large
- double distance_max = max_overlap_rate * rmin; // FIXME should be imposed by the user? for example max_gap_rate
+ double distance_max = max_overlap_rate * rmin;
if ( ( distance1 > distance_max)
|| ( distance2 > distance_max)
|| ( distance3 > distance_max)
@@ -577,7 +601,7 @@
void SpherePadder::place_at_tetra_centers ()
{
- cerr << "place at tetra centers... ";
+ BEGIN_FUNCTION("Place at tetra centers");
Sphere S;
S.type = AT_TETRA_CENTER;
@@ -610,12 +634,12 @@
place_sphere_4contacts(n);
}
- cerr << "Done" << endl;
+ END_FUNCTION;
}
void SpherePadder::place_at_tetra_vertexes ()
{
- cerr << "place at tetra vertexes... ";
+ BEGIN_FUNCTION("Place at tetra vertexes");
Sphere S;
S.type = AT_TETRA_VERTEX;
@@ -639,7 +663,7 @@
S.R = rmin;
- double pondere = 1.0 / 3.0; // FIXME parametrable
+ double pondere = .333333333; // FIXME parametrable
for (unsigned int n = 0 ; n < 4 ; ++n)
{
S.x = pondere * mesh->node[ T.nodeId[n] ].x + (1.0-pondere) * centre[0];
@@ -657,7 +681,7 @@
place_sphere_4contacts(n);
}
- cerr << "Done" << endl;
+ END_FUNCTION;
}
Modified: trunk/extra/SpherePadder/SpherePadder.hpp
===================================================================
--- trunk/extra/SpherePadder/SpherePadder.hpp 2009-01-09 13:09:26 UTC (rev 1618)
+++ trunk/extra/SpherePadder/SpherePadder.hpp 2009-01-09 14:44:43 UTC (rev 1619)
@@ -13,6 +13,9 @@
#include "TetraMesh.hpp"
+# define BEGIN_FUNCTION(arg) if (trace_functions) cerr << (arg) << "... "
+# define END_FUNCTION if (trace_functions) cerr << "Done\n"
+
enum SphereType {AT_NODE, AT_SEGMENT, AT_FACE, AT_TETRA_CENTER, AT_TETRA_VERTEX};
struct Sphere
@@ -35,41 +38,46 @@
class SpherePadder
{
-protected:
-
- vector<vector<unsigned int> > combination;
+ protected:
+
+ vector<vector<unsigned int> > combination;
- double distance_spheres (unsigned int i, unsigned int j);
- double distance_centre_spheres(Sphere& S1, Sphere& S2);
- double distance_vector3 (double V1[],double V2[]);
- void place_at_nodes ();
- void place_at_segment_middle ();
- void place_at_faces ();
- void place_at_tetra_centers ();
- void place_at_tetra_vertexes ();
- void cancel_overlap ();
- unsigned int place_fifth_sphere(unsigned int s1, unsigned int s2, unsigned int s3, unsigned int s4, Sphere& S);
- unsigned int place_sphere_4contacts (unsigned int sphereId);
-
- double rmin,rmax,rmoy,dr;
- double ratio;
- double max_overlap_rate;
- unsigned int n1,n2,n3,n4,n5,n_densify;
-
- TetraMesh * mesh;
- vector<Sphere> sphere;
-
- public:
+ double distance_spheres (unsigned int i, unsigned int j);
+ double distance_centre_spheres(Sphere& S1, Sphere& S2);
+ double distance_vector3 (double V1[],double V2[]);
+ void place_at_nodes ();
+ void place_at_segment_middle ();
+ void place_at_faces ();
+ void place_at_tetra_centers ();
+ void place_at_tetra_vertexes ();
+ void cancel_overlaps ();
+ unsigned int place_fifth_sphere(unsigned int s1, unsigned int s2, unsigned int s3, unsigned int s4, Sphere& S);
+ unsigned int place_sphere_4contacts (unsigned int sphereId);
+
+ double rmin,rmax,rmoy,dr;
+ double ratio;
+ double max_overlap_rate;
+ unsigned int n1,n2,n3,n4,n5,n_densify;
+ unsigned int nb_iter_max;
+
+ TetraMesh * mesh;
+ vector <Sphere> sphere;
+
+ bool trace_functions;
+
+ public:
- void plugTetraMesh (TetraMesh * mesh);
- void save_mgpost (const char* name);
- void save_Rxyz (const char* name);
-
- SpherePadder();
- // TODO destructor that clean TetraMesh*
+ bool meshIsPlugged;
+
+ void plugTetraMesh (TetraMesh * mesh);
+ void save_mgpost (const char* name);
+ void save_Rxyz (const char* name);
- void pad_5 ();
- // void densify ();
+ SpherePadder();
+ // TODO destructor that clean TetraMesh*?
+
+ void pad_5 ();
+ // void densify ();
};
Modified: trunk/extra/SpherePadder/TetraMesh.cpp
===================================================================
--- trunk/extra/SpherePadder/TetraMesh.cpp 2009-01-09 13:09:26 UTC (rev 1618)
+++ trunk/extra/SpherePadder/TetraMesh.cpp 2009-01-09 14:44:43 UTC (rev 1619)
@@ -10,8 +10,93 @@
#include "TetraMesh.hpp"
-void TetraMesh::read_data (const char* name)
+TetraMesh::TetraMesh ()
{
+ 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")
+ {
+ 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;
+ node.push_back(N);
+ }
+ }
+
+ if (token == "$Elements")
+ {
+ unsigned int nbElements;
+ unsigned int num_element, element_type, nbTags ;
+ Tetraedre T;
+ unsigned int t = 0;
+
+ meshFile >> nbElements;
+ for (unsigned int e = 0 ; e < nbElements ; ++e)
+ {
+ meshFile >> num_element >> element_type;
+ if (element_type != 4) // 4-node tetrahedron
+ {
+ meshFile.getline(not_read,150);
+ continue;
+ }
+
+
+ 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 >> T.nodeId[0] >> T.nodeId[1] >> T.nodeId[2] >> T.nodeId[3];
+
+ // 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;
+
+ 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);
+
+ tetraedre.push_back(T);
+ ++t;
+ }
+ }
+
+ if (token == "$EndElements") break;
+
+ meshFile >> token;
+ }
+
+ organize ();
+}
+
+void TetraMesh::read (const char* name)
+{
ifstream meshFile(name);
if(!meshFile)
{
@@ -19,7 +104,6 @@
return;
}
- cout << "Read data... " << flush;
string token;
meshFile >> token;
@@ -49,7 +133,7 @@
meshFile >> T.nodeId[0] >> T.nodeId[1] >> T.nodeId[2] >> T.nodeId[3];
// numbers begin at 0 instead of 1
- // (0 in C corresponds to 1 in the file)
+ // (0 in C/C++ corresponds to 1 in the file)
T.nodeId[0] -= 1;
T.nodeId[1] -= 1;
T.nodeId[2] -= 1;
@@ -68,7 +152,6 @@
meshFile >> token;
}
- cout << "Done" << endl;
organize ();
}
@@ -80,7 +163,7 @@
void TetraMesh::organize ()
{
- cout << "Organize data... " << flush;
+ //cout << "Organize data... " << flush;
// Translate all nodes in such a manner that all coordinates are > 0
xtrans = node[0].x;
@@ -92,9 +175,9 @@
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;
+ 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;
@@ -109,7 +192,7 @@
}
// Face creation
- vector <Face> tmpFace; // This will contain all faces more than one time
+ vector <Face> tmpFace; // This will contain all faces more than one time (with duplications)
Face F;
F.tetraOwner.push_back(0);
F.belongBoundary = true;
@@ -175,7 +258,7 @@
vector <Segment> tmpSegment;
Segment S;
S.faceOwner.push_back(0);
- //S.belongBoundary = true; // a voir
+
for (unsigned int i = 0 ; i < face.size() ; ++i)
{
S.faceOwner[0] = i;
@@ -251,7 +334,7 @@
if ( (tetraedre[t1].nodeId[0] > tetraedre[t2].nodeId[3])
|| (tetraedre[t1].nodeId[3] < tetraedre[t2].nodeId[0]) ) continue;
- // TODO mettre du while...
+ // FIXME mettre du while... (?)
for (unsigned int i = 0 ; i < 4 ; i++)
{
for (unsigned int j = 0 ; j < 4 ; j++)
@@ -267,17 +350,7 @@
}
}
- // Define Owners FIXME : a voir plus tard
- /*
- for (unsigned int t = 0 ; t < tetraedre.size() ; ++t)
- {
- node[tetraedre[t].nodeId[0]].tetraOwner.push_back(t);
- node[tetraedre[t].nodeId[1]].tetraOwner.push_back(t);
- node[tetraedre[t].nodeId[2]].tetraOwner.push_back(t);
- node[tetraedre[t].nodeId[3]].tetraOwner.push_back(t);
- }
- */
- cout << "Done" << endl;
+ isOrganized = true;
}
Modified: trunk/extra/SpherePadder/TetraMesh.hpp
===================================================================
--- trunk/extra/SpherePadder/TetraMesh.hpp 2009-01-09 13:09:26 UTC (rev 1618)
+++ trunk/extra/SpherePadder/TetraMesh.hpp 2009-01-09 14:44:43 UTC (rev 1619)
@@ -14,7 +14,7 @@
#include <vector>
#include <iostream>
#include <fstream>
-#include <stdlib.h> // for qsort
+#include <stdlib.h>
#include <math.h>
using namespace std;
@@ -54,22 +54,27 @@
{
protected:
- double xtrans,ytrans,ztrans;
- void organize ();
+
+ void organize ();
public:
- vector<Node> node;
- vector<Segment> segment;
- vector<Face> face;
- vector<Tetraedre> tetraedre;
+ vector <Node> node;
+ vector <Segment> segment;
+ vector <Face> face;
+ vector <Tetraedre> tetraedre;
double mean_segment_length;
double min_segment_length;
double max_segment_length;
- void read_data (const char* name);
-
+ TetraMesh();
+
+ double xtrans,ytrans,ztrans;
+ bool isOrganized;
+
+ void read (const char* name);
+ void read_gmsh (const char* name);
};
#endif // TETRA_MESH_HPP
Modified: trunk/extra/SpherePadder/main.cpp
===================================================================
--- trunk/extra/SpherePadder/main.cpp 2009-01-09 13:09:26 UTC (rev 1618)
+++ trunk/extra/SpherePadder/main.cpp 2009-01-09 14:44:43 UTC (rev 1619)
@@ -10,21 +10,231 @@
#include "SpherePadder.hpp"
+//#define DEBUG
+
+unsigned int mesh_format;
+vector <unsigned int> output_format;
+char mesh_file_name[100];
+char output_file_name[100];
+void resume();
+
int main()
-{
- TetraMesh * mesh = new TetraMesh();
- mesh->read_data("test.msh");
+{
+#ifdef DEBUG
+
+ TetraMesh * mesh = new TetraMesh();
+ mesh->read_gmsh("test2.msh");
- SpherePadder * padder = new SpherePadder();
- padder->plugTetraMesh(mesh);
-
- padder->pad_5();
-
- padder->save_mgpost("mgp.out.001");
- padder->save_Rxyz("out.Rxyz");
+ SpherePadder * padder = new SpherePadder();
+ padder->plugTetraMesh(mesh);
- return 0;
+ padder->pad_5();
+
+ padder->save_mgpost("mgp.out.001");
+ padder->save_Rxyz("out.Rxyz");
+
+ return 0;
+
+#else
+
+ char name_with_ext[100];
+
+ bool mesh_format_is_defined = false;
+ bool mesh_file_name_is_defined = false;
+ bool output_file_name_is_defined = false;
+
+ // default parameters
+ mesh_format = 1;
+ strcpy(mesh_file_name,"unknown");
+ strcpy(output_file_name,"unknown");
+
+ TetraMesh * mesh = new TetraMesh();
+ SpherePadder * padder = new SpherePadder();
+
+ unsigned int answer = 99;
+
+ while(1)
+ {
+ answer = 99;
+ cout << endl;
+ cout << "1. Define mesh format" << endl;
+ cout << "2. Define mesh file name" << endl;
+ cout << "3. Define output format" << endl;
+ cout << "4. Define output file name" << endl;
+ cout << "5. Resume" << endl;
+ cout << "6. GO" << endl;
+ cout << endl;
+ cout << "0. Quit" << endl;
+ cout << endl << "> ";
+ cin >> answer;
+
+ switch (answer)
+ {
+ case 0: return 0;
+
+ case 1:
+ answer = 99;
+ cout << endl;
+ cout << "1. SpherePadder" << endl;
+ cout << "2. gmsh" << endl;
+ cout << endl << "> ";
+ cin >> answer;
+ if (answer >= 1 && answer <=2)
+ {
+ mesh_format = answer;
+ mesh_format_is_defined = true;
+ }
+ break;
+
+ case 2:
+ cout << "Enter mesh file name: ";
+ cin >> mesh_file_name;
+ mesh_file_name_is_defined = true;
+ resume();
+ break;
+
+ case 3:
+ answer = 99;
+ cout << endl;
+ cout << "1. mgpost" << endl;
+ cout << "2. Rxyz" << endl;
+ cout << endl;
+ cout << "0. Done" << endl;
+ cout << endl << "> ";
+ cin >> answer;
+
+ if (answer == 0) break;
+ else if (answer >= 1 && answer <= 2)
+ {
+ bool added = false;
+ for (unsigned int i = 0 ; i < output_format.size() ; i++)
+ {
+ if (answer == output_format[i]) { added = true; break;}
+ }
+ if (!added) output_format.push_back(answer);
+ }
+ break;
+
+ case 4:
+ cout << "Enter output file name (without extension): ";
+ cin >> output_file_name;
+ output_file_name_is_defined = true;
+ resume();
+ break;
+
+ case 5: resume(); break;
+
+ case 6:
+
+ if (!mesh_file_name_is_defined)
+ {
+ cout << "mesh file name is not defined!" << endl;
+ break;
+ }
+
+ if (!output_file_name_is_defined)
+ {
+ cout << "output file name is not defined!" << endl;
+ break;
+ }
+
+ if (mesh_format_is_defined)
+ {
+ switch (mesh_format)
+ {
+ case 1:
+ mesh->read(mesh_file_name);
+ break;
+
+ case 2:
+ mesh->read_gmsh(mesh_file_name);
+ break;
+ }
+ }
+ else
+ {
+ cout << "mesh format is not defined!" << endl;
+ break;
+ }
+
+ if (!(padder->meshIsPlugged)) padder->plugTetraMesh(mesh);
+
+ if (output_format.empty())
+ {
+ cout << "output format not defined!" << endl;
+ break;
+ }
+
+ padder->pad_5(); // TODO controlled by user...
+
+ for (unsigned int i = 0 ; i < output_format.size() ; i++)
+ {
+ switch (output_format[i])
+ {
+ case 1:
+ strcpy(name_with_ext,output_file_name);
+ strcat(name_with_ext,".mgp");
+ padder->save_mgpost(name_with_ext);
+ cout << "file " << name_with_ext << " has been created" << endl;
+ break;
+
+ case 2:
+ strcpy(name_with_ext,output_file_name);
+ strcat(name_with_ext,".Rxyz");
+ padder->save_Rxyz(name_with_ext);
+ cout << "file " << name_with_ext << " has been created" << endl;
+ break;
+ }
+ }
+
+ return 0;
+
+ default:
+ break;
+
+ }
+ }
+
+ return 0;
+
+#endif // DEBUG
}
+void resume()
+{
+ cout << "--------------------------------------------------------" << endl;
+ cout << "mesh format: ";
+ switch (mesh_format)
+ {
+ case 1: cout << "SpherePadder" << endl; break;
+ case 2: cout << "gmsh" << endl;break;
+ default: cout << "unknown" << endl;
+ }
+
+ cout << "mesh file name: " << mesh_file_name << endl;
+ cout << "output format: ";
+ if ((!output_format.empty()))
+ {
+ cout << "[";
+ for (unsigned int i = 0 ; i < output_format.size() ; i++)
+ {
+ switch (output_format[i])
+ {
+ case 1: cout << " mgpost ";break;
+ case 2: cout << " Rxyz ";break;
+ default: break;
+ }
+ }
+ cout << "]" << endl;
+ }
+ else cout << "unknown" << endl;
+
+ cout << "output file name (without extension): " << output_file_name << endl;
+
+ cout << "--------------------------------------------------------" << endl;
+}
+
+
+
