yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #01187
[svn] r1747 - in trunk/extra/SpherePadder: . SpherePackingTriangulation
Author: richefeu
Date: 2009-04-06 17:42:25 +0200 (Mon, 06 Apr 2009)
New Revision: 1747
Modified:
trunk/extra/SpherePadder/CellPartition.cpp
trunk/extra/SpherePadder/CellPartition.hpp
trunk/extra/SpherePadder/SpherePackingTriangulation/SpherePackingTriangulation.cpp
trunk/extra/SpherePadder/SpherePackingTriangulation/SpherePackingTriangulation.hpp
trunk/extra/SpherePadder/SpherePadder.cpp
trunk/extra/SpherePadder/SpherePadder.hpp
trunk/extra/SpherePadder/main.cpp
Log:
- overlap with spheres inserted by user -> bug fixed
- Densification begins to work
Modified: trunk/extra/SpherePadder/CellPartition.cpp
===================================================================
--- trunk/extra/SpherePadder/CellPartition.cpp 2009-04-06 10:26:50 UTC (rev 1746)
+++ trunk/extra/SpherePadder/CellPartition.cpp 2009-04-06 15:42:25 UTC (rev 1747)
@@ -45,7 +45,7 @@
ksize = (unsigned int)((zmax - zmin) / (security_factor * mesh.mean_segment_length));
if (ksize < 1) ksize = 1;
- //cerr << "nb cells: " << isize << ", " << jsize << ", " << ksize << endl;
+ cout << "nb cells: " << isize << ", " << jsize << ", " << ksize << endl;
vector<unsigned int> kvec;
for (unsigned int k = 0 ; k < ksize ; ++k) kvec.push_back(0);
@@ -93,8 +93,12 @@
cell[ cellId[current_i][current_j][current_k] ].sphereId.push_back(n);
}
-
+void CellPartition::add_in_cell(unsigned int n, unsigned int i, unsigned int j, unsigned int k)
+{
+ cell[ cellId[i][j][k] ].sphereId.push_back(n);
+}
+
void CellPartition::locateCellOf(double x, double y, double z)
{
int i,j,k;
Modified: trunk/extra/SpherePadder/CellPartition.hpp
===================================================================
--- trunk/extra/SpherePadder/CellPartition.hpp 2009-04-06 10:26:50 UTC (rev 1746)
+++ trunk/extra/SpherePadder/CellPartition.hpp 2009-04-06 15:42:25 UTC (rev 1747)
@@ -39,6 +39,7 @@
CellPartition();
void init(TetraMesh & mesh, double security_factor = 1.0);
void add(unsigned int n, double x, double y, double z);
+ void add_in_cell(unsigned int n, unsigned int i, unsigned int j, unsigned int k);
void locateCellOf(double x, double y, double z);
Cell& get_cell (unsigned int i,unsigned int j,unsigned int k) { return cell[ cellId[i][j][k] ]; }
Modified: trunk/extra/SpherePadder/SpherePackingTriangulation/SpherePackingTriangulation.cpp
===================================================================
--- trunk/extra/SpherePadder/SpherePackingTriangulation/SpherePackingTriangulation.cpp 2009-04-06 10:26:50 UTC (rev 1746)
+++ trunk/extra/SpherePadder/SpherePackingTriangulation/SpherePackingTriangulation.cpp 2009-04-06 15:42:25 UTC (rev 1747)
@@ -61,15 +61,43 @@
return true;
}
-
-
float SpherePackingTriangulation::current_tetrahedron_get_volume()
{
Cell_handle cell = tetrahedron_iterator;
return (Tetrahedron(cell->vertex(0)->point(), cell->vertex(1)->point(),
cell->vertex(2)->point(), cell->vertex (3)->point())).volume();
}
-
+
+// warning order of radii is the same as id
+void SpherePackingTriangulation::current_tetrahedron_get_circumcenter(Real R1,Real R2,Real R3,Real R4,Real& x, Real& y, Real& z)
+{
+ Cell_handle cell = tetrahedron_iterator;
+
+ Real p1x = cell->vertex(0)->point().x();
+ Real p1y = cell->vertex(0)->point().y();
+ Real p1z = cell->vertex(0)->point().z();
+
+ Real p2x = cell->vertex(1)->point().x();
+ Real p2y = cell->vertex(1)->point().y();
+ Real p2z = cell->vertex(1)->point().z();
+
+ Real p3x = cell->vertex(2)->point().x();
+ Real p3y = cell->vertex(2)->point().y();
+ Real p3z = cell->vertex(2)->point().z();
+
+ Real p4x = cell->vertex(3)->point().x();
+ Real p4y = cell->vertex(3)->point().y();
+ Real p4z = cell->vertex(3)->point().z();
+
+ CGAL::weighted_circumcenterC3 (
+ p1x, p1y, p1z, R1*R1,
+ p2x, p2y, p2z, R2*R2,
+ p3x, p3y, p3z, R3*R3,
+ p4x, p4y, p4z, R4*R4,
+ x, y, z
+ );
+}
+
void SpherePackingTriangulation::current_tetrahedron_get_nodes(unsigned int & id1, unsigned int & id2, unsigned int & id3, unsigned int & id4)
{
Cell_handle cell = tetrahedron_iterator;
Modified: trunk/extra/SpherePadder/SpherePackingTriangulation/SpherePackingTriangulation.hpp
===================================================================
--- trunk/extra/SpherePadder/SpherePackingTriangulation/SpherePackingTriangulation.hpp 2009-04-06 10:26:50 UTC (rev 1746)
+++ trunk/extra/SpherePadder/SpherePackingTriangulation/SpherePackingTriangulation.hpp 2009-04-06 15:42:25 UTC (rev 1747)
@@ -83,6 +83,13 @@
public:
SpherePackingTriangulation();
+ bool has_no_cells() { return (tri.number_of_finite_cells() == 0); }
+ void clear()
+ {
+ tri.clear();
+ volumeAreComputed = false;
+ id_vh_link.clear();
+ }
bool insert_node(Real x, Real y, Real z, unsigned int id, bool isVirtual = false);
@@ -94,7 +101,28 @@
float current_tetrahedron_get_volume();
void current_tetrahedron_get_nodes(unsigned int & id1, unsigned int & id2, unsigned int & id3, unsigned int & id4);
-
+ void current_tetrahedron_get_circumcenter(Real R1,Real R2,Real R3,Real R4,Real& x, Real& y, Real& z);
+
+// TODO
+/*
+void current_tetrahedron_get_weighted_circumcenter(R1,R2,R3,R4,&x,&y,&z)
+{
+ const CGAL_Sphere& S0 = cell->vertex(0)->point();
+ const CGAL_Sphere& S1 = cell->vertex(1)->point();
+ const CGAL_Sphere& S2 = cell->vertex(2)->point();
+ const CGAL_Sphere& S3 = cell->vertex(3)->point();
+ Real x,y,z;
+
+ CGAL::weighted_circumcenterC3 (
+ S0.point().x(), S0.point().y(), S0.point().z(), S0.weight(),
+ S1.point().x(), S1.point().y(), S1.point().z(), S1.weight(),
+ S2.point().x(), S2.point().y(), S2.point().z(), S2.weight(),
+ S3.point().x(), S3.point().y(), S3.point().z(), S3.weight(),
+ x, y, z
+ );
+}
+*/
+
};
Modified: trunk/extra/SpherePadder/SpherePadder.cpp
===================================================================
--- trunk/extra/SpherePadder/SpherePadder.cpp 2009-04-06 10:26:50 UTC (rev 1746)
+++ trunk/extra/SpherePadder/SpherePadder.cpp 2009-04-06 15:42:25 UTC (rev 1747)
@@ -21,7 +21,7 @@
{
double d1 = (*(tetra_porosity *)a).void_volume;
double d2 = (*(tetra_porosity *)b).void_volume;
- return (d1 < d2) ? 1 :-1;
+ return ((d1) < (d2)) ? 1 :-1;
}
int compareDouble (const void * a, const void * b)
@@ -33,126 +33,106 @@
{
vector <unsigned int> lst;
unsigned int nb = 5;
-
+
for (unsigned int i = 0 ; i <= nb ; ++i)
- for (unsigned int j = i+1 ; j <= nb+1 ; ++j)
- for (unsigned int k = j+1 ; k <= nb+2 ; ++k)
- for (unsigned int l = k+1 ; l <= nb+3 ; ++l)
{
- lst.clear();
- lst.push_back(i);
- lst.push_back(j);
- lst.push_back(k);
- lst.push_back(l);
- combination.push_back(lst);
+ for (unsigned int j = i+1 ; j <= nb+1 ; ++j)
+ {
+ for (unsigned int k = j+1 ; k <= nb+2 ; ++k)
+ {
+ for (unsigned int l = k+1 ; l <= nb+3 ; ++l)
+ {
+ lst.clear();
+ lst.push_back(i);
+ lst.push_back(j);
+ lst.push_back(k);
+ lst.push_back(l);
+ combination.push_back(lst);
+ }
+ }
+ }
}
-
- max_overlap_rate = 1e-4;
- n1 = n2 = n3 = n4 = n5 = n_densify = 0;
- trace_functions = true;
- meshIsPlugged = false;
- probeIsDefined = false;
- ratio = 4.0; rmoy = 0.0;
- virtual_radius_factor = 5.0;
-
-/* FIXME
- pour le moment, l'utilisateur ne peut entre qu'un ratio.
- Le rayon des sphere depend du maillage (des longueurs des segments)
-*/
-
-}
-
-void SpherePadder::check_inProbe(unsigned int i)
-{
- if (!probeIsDefined) return;
- if (sphere[i].R <= 0.0) return;
- double dx = sphere[i].x - xProbe;
- double dy = sphere[i].y - yProbe;
- double dz = sphere[i].z - zProbe;
- double squared_distance = dx*dx + dy*dy + dz*dz;
- if (squared_distance + sphere[i].R*sphere[i].R < RProbe * RProbe) sphereInProbe.push_back(i);
+ n1 = n2 = n3 = n4 = n5 = n_densify = 0;
+ trace_functions = true;
+ meshIsPlugged = false;
+ probeIsDefined = false;
+ RadiusDataIsOK = RadiusIsSet = false;
+
+ max_overlap_rate = 1e-4;
+ virtual_radius_factor = 100.0;
}
-void SpherePadder::add_spherical_probe(double Rfact)
+void SpherePadder::setRadiusRatio(double r)
{
- double xmin,xmax;
- double ymin,ymax;
- double zmin,zmax;
- double R;
+ r = fabs(r);
+ if (r < 1.0) ratio = 1.0/r;
+ else ratio = r;
- xmin=xmax=mesh->node[0].x;
- ymin=ymax=mesh->node[0].y;
- zmin=zmax=mesh->node[0].z;
- for (unsigned int i = 1 ; i < mesh->node.size() ; ++i)
+ if (meshIsPlugged)
{
- xmin = (xmin > mesh->node[i].x) ? mesh->node[i].x : xmin;
- xmax = (xmax < mesh->node[i].x) ? mesh->node[i].x : xmax;
- ymin = (ymin > mesh->node[i].y) ? mesh->node[i].y : ymin;
- ymax = (ymax < mesh->node[i].y) ? mesh->node[i].y : ymax;
- zmin = (zmin > mesh->node[i].z) ? mesh->node[i].z : zmin;
- zmax = (zmax < mesh->node[i].z) ? mesh->node[i].z : zmax;
+ rmoy = 0.125 * mesh->mean_segment_length; // 1/8 = 0.125
+ rmin = (2.0 * rmoy) / (ratio + 1.0);
+ rmax = 2.0 * rmoy - rmin;
+ dr = rmax - rmoy;
+ RadiusDataIsOK = true;
+ cout << "rmin = " << rmin << endl;
+ cout << "rmax = " << rmax << endl;
}
- xProbe = 0.5 * (xmin + xmax);
- yProbe = 0.5 * (ymin + ymax);
- zProbe = 0.5 * (zmin + zmax);
- RProbe = (xmax - xmin);
- RProbe = (RProbe < (R = ymax - ymin)) ? RProbe : R;
- RProbe = (RProbe < (R = zmax - zmin)) ? RProbe : R;
- RProbe *= 0.5 * Rfact;
- probeIsDefined = true;
- compacity_file.open("compacity.dat");
+ else
+ {
+ rmin = rmax = rmoy = 0.0;
+ RadiusDataIsOK = false;
+ // rmin, rmax and rmoy will be calculated in plugTetraMesh
+ }
+ RadiusIsSet = true;
}
-
-double SpherePadder::compacity_in_probe(unsigned int ninsered)
+
+void SpherePadder::setRadiusRange(double min, double max)
{
- if (!probeIsDefined) return -1.0;
-
- double dr = 0.0;
- double Vs = 0.0;
- double fact = 1.333333333 * M_PI;
- for (unsigned int i = 1 ; i < sphereInProbe.size() ; ++i)
+ if (min > max)
{
- Vs += fact * sphere[i].R * sphere[i].R * sphere[i].R;
+ rmin = max;
+ rmax = min;
}
-
- // TODO ameliorer cette fonction en prennant en compte les spheres coupees...
-
- dr = Vs / (fact*RProbe*RProbe*RProbe);
- compacity_file << ninsered << ' ' << dr << endl;
- return (dr);
+ else
+ {
+ rmin = min;
+ rmax = max;
+ }
+ ratio = rmax/rmin;
+ rmoy = 0.5*(rmin+rmax);
+ RadiusDataIsOK = true;
+ RadiusIsSet = true;
}
+
void SpherePadder::plugTetraMesh (TetraMesh * pluggedMesh)
{
mesh = pluggedMesh;
partition.init(*mesh);
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 = 0.125
- rmin = (2.0 * rmoy) / (ratio + 1.0);
- rmax = 2.0 * rmoy - rmin;
- dr = rmax - rmoy;
- cerr << "rmax = " << rmax << endl;
- }
+
+ cout << "mesh->mean_segment_length = " << mesh->mean_segment_length << endl;
+ cout << "mesh->min_segment_length = " << mesh->min_segment_length << endl;
+ cout << "mesh->max_segment_length = " << mesh->max_segment_length << endl;
+
+ if (!RadiusDataIsOK && RadiusIsSet && ratio != 0.0) setRadiusRatio(ratio);
+ else cerr << "@SpherePadder::plugTetraMesh, no radius range defined!" << endl;
+
}
void SpherePadder::pad_5 ()
{
if (mesh == 0)
{
- cerr << "SpherePadder::pad_5, no mesh defined!" << endl;
+ cerr << "@SpherePadder::pad_5, no mesh defined!" << endl;
return;
}
if (!(mesh->isOrganized))
{
- cerr << "SpherePadder::pad_5, mesh is not valid!" << endl;
+ cerr << "@SpherePadder::pad_5, mesh is not valid!" << endl;
return;
}
@@ -187,86 +167,229 @@
}
-void SpherePadder::densify() // TODO FIXME loop until target solid fraction is reached
+void SpherePadder::densify()
{
BEGIN_FUNCTION ("Densify");
+ unsigned int added = 0, back_added = 0;
+ unsigned int nbfail = 0;
+ unsigned int max_iter_densify = 40;
+ unsigned int i;
+
+ repack_null_radii();
+
+ for (i = 0 ; i < max_iter_densify ; ++i)
+ {
+ // CRITERE BIDON !!!!! WARNING
+ if ((added = iter_densify()) == back_added)
+ {
+ if (++nbfail == 3) { cout << "@densify, cannot add more spheres with this ratio" << endl; break; }
+ }
+ else nbfail = 0;
+
+ cout << "iter " << i << ", nb spheres = " << sphere.size() /*- nzero*/ << ", added = " << added << endl;
+ back_added = added;
+ }
+
+ if (i == max_iter_densify) cout << "@densify, maximum number of iteration reached" << endl;
+
+ END_FUNCTION;
+}
+
+void SpherePadder::repack_null_radii() // TODO FIXME loop until target solid fraction is reached
+{
+ for (unsigned int i = 0 ; i < sphere.size() ; i++)
+ {
+ if (sphere[i].R > 0.0) continue;
+ place_sphere_4contacts(i);
+ }
+}
+
+unsigned int SpherePadder::iter_densify(unsigned int nb_check)
+{
+ unsigned int nb_added = 0, total_added = 0;
tetra_porosity P;
- Sphere S,S1,S2,S3,S4;
- bool can_be_added;
- S.type = INSERTED_BY_USER;//FROM_TRIANGULATION;
- unsigned int ns = sphere.size();
+ Sphere S;
+ S.type = FROM_TRIANGULATION;
+ unsigned int ns = sphere.size();
+ unsigned int failure;
+ unsigned int nbVirtual;
+ //double Volume_min = 1.333333333*M_PI*rmin*rmin*rmin;
- //place_virtual_spheres();
+ // TODO add bool VirtualSpheresOK
- // TODO clear triangulation
+ triangulation.clear();
for (unsigned int i = 0 ; i < sphere.size() ; i++)
{
if (sphere[i].R <= 0.0) continue;
- triangulation.insert_node(sphere[i].x, sphere[i].y, sphere[i].z, i, false);
+ triangulation.insert_node(sphere[i].x, sphere[i].y, sphere[i].z, i, (sphere[i].type == VIRTUAL));
}
-
- if (!tetra_porosities.empty()) tetra_porosities.clear();
- // FIXME if triangulation is empty ...
+ if (triangulation.has_no_cells()) return total_added;
+
+ tetra_porosities.clear();
triangulation.init_current_tetrahedron();
do
{
triangulation.current_tetrahedron_get_nodes(P.id1,P.id2,P.id3,P.id4);
- if (P.id1>=ns || P.id2>=ns || P.id3>=ns || P.id4>=ns) continue;
+ if (P.id1 >= ns || P.id2 >= ns || P.id3 >= ns || P.id4 >= ns) continue; // FIXME why ?
+
+ nbVirtual = 0;
+ if (sphere[P.id1].type == VIRTUAL) ++nbVirtual;
+ if (sphere[P.id2].type == VIRTUAL) ++nbVirtual;
+ if (sphere[P.id3].type == VIRTUAL) ++nbVirtual;
+ if (sphere[P.id4].type == VIRTUAL) ++nbVirtual;
+ if (nbVirtual == 4) continue;
+
+/*
+ if (squared_distance_centre_spheres(P.id1, P.id2) > 0.01) continue;
+ if (squared_distance_centre_spheres(P.id1, P.id3) > 0.01) continue;
+ if (squared_distance_centre_spheres(P.id1, P.id4) > 0.01) continue;
+ if (squared_distance_centre_spheres(P.id2, P.id3) > 0.01) continue;
+ if (squared_distance_centre_spheres(P.id2, P.id4) > 0.01) continue;
+ if (squared_distance_centre_spheres(P.id3, P.id4) > 0.01) continue;
+*/
+
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]);
- if (P.void_volume>0.0)
+ if (P.void_volume > /*Volume_min*/ 0.0 || nbVirtual > 0)
{
tetra_porosities.push_back(P);
}
} while (triangulation.next_tetrahedron());
- // sort tetrahdrons from big to small void volumes
+ // sort tetrahdrons from bigger to smaller void volumes
qsort(&(tetra_porosities[0]),tetra_porosities.size(),sizeof(tetra_porosity),compare_tetra_porosity);
-
- // TODO FIXME exclude tetra that involve virtual spheres
for (unsigned int i = 0 ; i < tetra_porosities.size() ; i++)
- {
- //cout << tetra_porosities[i].volume << "\t" << tetra_porosities[i].void_volume << endl;
- //if (tetra_porosities[i].void_volume > 0.0)
+ {
+ failure = place_fifth_sphere(tetra_porosities[i].id1,tetra_porosities[i].id2,tetra_porosities[i].id3,tetra_porosities[i].id4,S);
- S1 = sphere[tetra_porosities[i].id1];
- S2 = sphere[tetra_porosities[i].id2];
- S3 = sphere[tetra_porosities[i].id3];
- S4 = sphere[tetra_porosities[i].id4];
+ if (!failure)
+ {
+ failure = check_overlaps(S,ns+1); // ns+1 i.e. on exclut aucune sphere
+ }
+ else // failure
+ {
+// Sphere & S1 = sphere[tetra_porosities[i].id1];
+// Sphere & S2 = sphere[tetra_porosities[i].id2];
+// Sphere & S3 = sphere[tetra_porosities[i].id3];
+// Sphere & S4 = sphere[tetra_porosities[i].id4];
- can_be_added = (S1.type != VIRTUAL);
- can_be_added &= (S2.type != VIRTUAL);
- can_be_added &= (S3.type != VIRTUAL);
- can_be_added &= (S4.type != VIRTUAL);
+ /*
+ CGAL::weighted_circumcenterC3 (
+ S1.x, S1.y, S1.z, S1.R * S1.R,
+ S2.x, S2.y, S2.z, S2.R * S2.R,
+ S3.x, S3.y, S3.z, S3.R * S3.R,
+ S4.x, S4.y, S4.z, S4.R * S4.R,
+ S.x, S.y, S.z
+ );
+ */
+
+
+
+// double distance_min = distance_centre_spheres(S,S1) - S1.R;
+// double distance;
+// distance_min = (distance_min < (distance = distance_centre_spheres(S,S2) - S2.R)) ? distance_min : distance;
+// distance_min = (distance_min < (distance = distance_centre_spheres(S,S3) - S3.R)) ? distance_min : distance;
+// distance_min = (distance_min < (distance = distance_centre_spheres(S,S4) - S4.R)) ? distance_min : distance;
+// S.R = distance_min;
+// if (S.R > rmax) S.R = rmax;
+ // TODO mettre en contact
+
+ S.R = rmin; // test
+
+ //if (S.R >= rmin)
+ {
+//failure = 0;
+ failure = check_overlaps(S,ns+10);
+ }
+ }
- if (can_be_added)
+ if(!failure)
{
-// TODO write a function 'try_to_place_into(id1,id2,id3,id4)'
+ vector<neighbor_with_distance> neighbor;
+ build_sorted_list_of_neighbors(S, neighbor);
- if( place_fifth_sphere(tetra_porosities[i].id1,tetra_porosities[i].id2,tetra_porosities[i].id3,tetra_porosities[i].id4,S) )
+ S.R += neighbor[0].distance;
+
+ if (S.R >= rmin && S.R <= rmax)
{
- sphere.push_back(S);
+ sphere.push_back(S);++(n_densify);
+ partition.add(ns,S.x,S.y,S.z);
+ triangulation.insert_node(S.x,S.y,S.z,ns,false); // necessary to compute mean solid fraction
+ ++ns; ++nb_added; ++total_added;
}
+ }
- /*
- S.x = 0.25 * (S1.x + S2.x + S3.x + S4.x);
- S.y = 0.25 * (S1.y + S2.y + S3.y + S4.y);
- S.z = 0.25 * (S1.z + S2.z + S3.z + S4.z);
- S.R = rmin;
-
- sphere.push_back(S);
- place_sphere_4contacts(ns,4);
- ++ns;
- */
+ // check if we must stop
+ if (nb_added >= nb_check)
+ {
+ //cout << "check, total_added = " << total_added << endl;
+ nb_added = 0;
+ // getMeanSolidFraction();
}
}
- END_FUNCTION;
+return total_added;
}
+void SpherePadder::save_tri_mgpost (const char* name)
+{
+ // triangulation
+ triangulation.clear();
+ for (unsigned int i = 0 ; i < sphere.size() ; i++)
+ {
+ if (sphere[i].R <= 0.0) continue;
+ triangulation.insert_node(sphere[i].x, sphere[i].y, sphere[i].z, i, (sphere[i].type == VIRTUAL));
+ }
+ if (triangulation.has_no_cells()) return;
+ ofstream fmgpost(name);
+ fmgpost << "<?xml version=\"1.0\"?>" << endl
+ << " <mgpost mode=\"3D\">" << endl
+ << " <state id=\"" << 1 << "\" time=\"" << 0.0 << "\">" << endl;
+
+ double xp,yp,zp,rad,void_volume;
+ unsigned int id1,id2,id3,id4,i=0;
+ tetra_porosities.clear();
+ triangulation.init_current_tetrahedron();
+ do
+ {
+ triangulation.current_tetrahedron_get_nodes(id1,id2,id3,id4);
+ //if(sphere[id1].type == VIRTUAL && sphere[id2].type == VIRTUAL && sphere[id3].type == VIRTUAL && sphere[id4].type == VIRTUAL) continue;
+
+ xp = 0.25*(sphere[id1].x+sphere[id2].x+sphere[id3].x+sphere[id4].x);
+ yp = 0.25*(sphere[id1].y+sphere[id2].y+sphere[id3].y+sphere[id4].y);
+ zp = 0.25*(sphere[id1].z+sphere[id2].z+sphere[id3].z+sphere[id4].z);
+ rad = (sphere[id1].R>sphere[id2].R)?sphere[id1].R:sphere[id2].R;
+ rad = (rad>sphere[id3].R)?rad:sphere[id3].R;
+ rad = (rad>sphere[id4].R)?rad:sphere[id4].R;
+
+ void_volume = triangulation.current_tetrahedron_get_volume() - solid_volume_of_tetrahedron(sphere[id1], sphere[id2], sphere[id3], sphere[id4]);
+ if (void_volume<0.0) void_volume = 0.0;
+
+ fmgpost << " <body>" << endl;
+ fmgpost << " <POLYE id=\"" << ++i << "\" r=\"" << rad << "\">" << endl
+ << " <position x=\"" << xp << "\" y=\"" << yp << "\" z=\"" << zp << "\"/>" << endl
+ //<< " <velocity rot=\"" << void_volume << "\"/>" << endl
+ << " <velocity rot=\"" << triangulation.current_tetrahedron_get_volume() << "\"/>" << endl
+ << " <node x=\""<< sphere[id1].x-xp << "\" y=\"" << sphere[id1].y-yp << "\" z=\"" << sphere[id1].z-zp << "\"/>" << endl
+ << " <node x=\""<< sphere[id2].x-xp << "\" y=\"" << sphere[id2].y-yp << "\" z=\"" << sphere[id2].z-zp << "\"/>" << endl
+ << " <node x=\""<< sphere[id3].x-xp << "\" y=\"" << sphere[id3].y-yp << "\" z=\"" << sphere[id3].z-zp << "\"/>" << endl
+ << " <node x=\""<< sphere[id4].x-xp << "\" y=\"" << sphere[id4].y-yp << "\" z=\"" << sphere[id4].z-zp << "\"/>" << endl
+ << " <face n1=\"0\" n2=\"1\" n3=\"2\"/>" << endl
+ << " <face n1=\"0\" n2=\"2\" n3=\"3\"/>" << endl
+ << " <face n1=\"0\" n2=\"3\" n3=\"1\"/>" << endl
+ << " <face n1=\"1\" n2=\"3\" n3=\"2\"/>" << endl
+ << " </POLYE>" << endl << flush;
+ fmgpost << " </body>" << endl;
+ } while (triangulation.next_tetrahedron());
+
+ fmgpost << " </state>" << endl
+ << " </mgpost>" << endl;
+
+}
+
void SpherePadder::save_mgpost (const char* name)
{
BEGIN_FUNCTION ("Save mgp");
@@ -279,19 +402,22 @@
fmgpost << "<?xml version=\"1.0\"?>" << endl
<< " <mgpost mode=\"3D\">" << endl
+ // Note that number of newcolors is limited in mgpost
<< " <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
<< " <newcolor name=\"insered by user\"/>" << endl
- //<< " <newcolor name=\"from triangulation\"/>" << endl
+ << " <newcolor name=\"from triangulation\"/>" << endl
<< " <state id=\"" << 1
<< "\" time=\"" << 0.0 << "\">" << endl;
for (unsigned int i = 0 ; i < sphere.size() ; ++i)
{
- if (sphere[i].R <= 0.0 && sphere[i].type != AT_NODE) continue;
+ // the sphere with R=0 that are placed at nodes are saved because mgpost needs them
+ // for the display of tetrahdrons.
+ if (sphere[i].R <= 0.0 && sphere[i].type != AT_NODE) continue;
if (sphere[i].type == VIRTUAL) continue;
@@ -424,7 +550,6 @@
if (S.R > rmax) S.R = rmoy + dr * (double)rand()/(double)RAND_MAX;
sphere.push_back(S); ++(n2);
- //check_inProbe(ns);compacity_in_probe(ns);
partition.add(ns,S.x,S.y,S.z);
mesh->segment[s].sphereId = ns;
@@ -471,16 +596,9 @@
sphere.push_back(S); ++(n3);
place_sphere_4contacts(ns);
- //partition.add(ns,S.x,S.y,S.z); // test
++ns;
}
-// for (unsigned int n = (n1+n2) ; n < sphere.size() ; ++n)
-// {
-// place_sphere_4contacts(n);
-// check_inProbe(n);compacity_in_probe(n);
-// }
-
END_FUNCTION;
}
@@ -510,19 +628,11 @@
S.R = rmin;
- //S.tetraOwner = t;
- //mesh->tetraedre[t].sphereId.push_back(ns++);
sphere.push_back(S); ++(n4);
place_sphere_4contacts(ns);
++ns;
}
-// for (unsigned int n = (n1+n2+n3) ; n < sphere.size() ; ++n)
-// {
-// place_sphere_4contacts(n);
-// check_inProbe(n);compacity_in_probe(n);
-// }
-
END_FUNCTION;
}
@@ -552,33 +662,32 @@
S.R = rmin;
- double pondere = 0.333333333; // FIXME parametrable
+ double pondere = 0.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];
S.y = pondere * mesh->node[ T.nodeId[n] ].y + (1.0 - pondere) * centre[1];
S.z = pondere * mesh->node[ T.nodeId[n] ].z + (1.0 - pondere) * centre[2];
- //S.tetraOwner = t;
- //mesh->tetraedre[t].sphereId.push_back(ns++);
sphere.push_back(S); ++(n5);
place_sphere_4contacts(ns);
++ns;
}
}
-// for (unsigned int n = (n1+n2+n3+n4) ; n < sphere.size() ; ++n)
-// {
-// place_sphere_4contacts(n);
-// check_inProbe(n);compacity_in_probe(n);
-// }
-//
END_FUNCTION;
}
-
+double SpherePadder::squared_distance_centre_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 ((lx*lx + ly*ly + lz*lz));
+}
-
+// deprecated
double SpherePadder::distance_spheres(unsigned int i, unsigned int j)
{
double lx,ly,lz;
@@ -588,6 +697,15 @@
return (sqrt(lx*lx + ly*ly + lz*lz) - sphere[i].R - sphere[j].R);
}
+double SpherePadder::distance_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) - S1.R - S2.R);
+}
+
double SpherePadder::distance_centre_spheres(Sphere& S1, Sphere& S2)
{
double lx,ly,lz;
@@ -597,44 +715,86 @@
return (sqrt(lx*lx + ly*ly + lz*lz));
}
-unsigned int SpherePadder::place_sphere_4contacts (unsigned int sphereId, unsigned int nb_combi_max)
-{
- Sphere S = sphere[sphereId];
- Sphere Sbackup;
-
- vector<neighbor_with_distance> neighbor;
+void SpherePadder::build_sorted_list_of_neighbors(Sphere & S, vector<neighbor_with_distance> & neighbor)
+{
neighbor_with_distance N;
-
+
unsigned int id;
Cell current_cell;
partition.locateCellOf(S.x,S.y,S.z);
+
+ // Build the local list of neighbors
+ for (unsigned int i = partition.i_down() ; i <= partition.i_up() ; ++i)
+ {
+ for (unsigned int j = partition.j_down() ; j <= partition.j_up() ; ++j)
+ {
+ for (unsigned int k = partition.k_down() ; k <= partition.k_up() ; ++k)
+ {
+
+ current_cell = partition.get_cell(i,j,k);
+ for (unsigned int s = 0 ; s < current_cell.sphereId.size() ; ++s)
+ {
+ id = current_cell.sphereId[s];
+ //if (id != sphereId)
+ {
+ N.sphereId = id;
+ N.distance = distance_spheres(S,sphere[id]);
+ neighbor.push_back(N);
+ }
+ }
+
+ }
+ }
+ }
+
+ qsort(&(neighbor[0]),neighbor.size(),sizeof(neighbor_with_distance),compare_neighbor_with_distance);
+}
- // Build the list of neighbors
+void SpherePadder::build_sorted_list_of_neighbors(unsigned int sphereId, vector<neighbor_with_distance> & neighbor)
+{
+ neighbor_with_distance N;
+
+ unsigned int id;
+ Cell current_cell;
+
+ partition.locateCellOf(sphere[sphereId].x,sphere[sphereId].y,sphere[sphereId].z);
+
+ // Build the local list of neighbors
for (unsigned int i = partition.i_down() ; i <= partition.i_up() ; ++i)
{
- for (unsigned int j = partition.j_down() ; j <= partition.j_up() ; ++j)
- {
- for (unsigned int k = partition.k_down() ; k <= partition.k_up() ; ++k)
- {
-
- current_cell = partition.get_cell(i,j,k);
- for (unsigned int s = 0 ; s < current_cell.sphereId.size() ; ++s)
- {
- id = current_cell.sphereId[s];
- if (id != sphereId)
- {
- N.sphereId = id;
- N.distance = distance_spheres(sphereId,id);
- neighbor.push_back(N);
- }
- }
-
- }
- }
+ for (unsigned int j = partition.j_down() ; j <= partition.j_up() ; ++j)
+ {
+ for (unsigned int k = partition.k_down() ; k <= partition.k_up() ; ++k)
+ {
+
+ current_cell = partition.get_cell(i,j,k);
+ for (unsigned int s = 0 ; s < current_cell.sphereId.size() ; ++s)
+ {
+ id = current_cell.sphereId[s];
+ if (id != sphereId)
+ {
+ N.sphereId = id;
+ N.distance = distance_spheres(sphereId,id);
+ neighbor.push_back(N);
+ }
+ }
+
+ }
+ }
}
-
- qsort(&(neighbor[0]),neighbor.size(),sizeof(neighbor_with_distance),compare_neighbor_with_distance);
+
+ qsort(&(neighbor[0]),neighbor.size(),sizeof(neighbor_with_distance),compare_neighbor_with_distance);
+}
+
+unsigned int SpherePadder::place_sphere_4contacts (unsigned int sphereId, unsigned int nb_combi_max)
+{
+ Sphere S = sphere[sphereId];
+ Sphere Sbackup;
+
+ vector<neighbor_with_distance> neighbor;
+ build_sorted_list_of_neighbors(sphereId, neighbor);
+
S.R += neighbor[0].distance;
if (S.R >= rmin && S.R <= rmax) sphere[sphereId].R = S.R;
else if (S.R > rmax) sphere[sphereId].R = rmax;
@@ -664,32 +824,10 @@
S = Sbackup;
failure = place_fifth_sphere(s1,s2,s3,s4,S);
-
+
if (!failure)
- {
- partition.locateCellOf(S.x,S.y,S.z);
- unsigned int nb = 0;
- for (unsigned int i = partition.i_down() ; i <= partition.i_up() ; ++i)
- {
- for (unsigned int j = partition.j_down() ; j <= partition.j_up() ; ++j)
- {
- for (unsigned int k = partition.k_down() ; k <= partition.k_up() ; ++k)
- {
- current_cell = partition.get_cell(i,j,k);
- nb += current_cell.sphereId.size(); // FIXME why ??
- for (unsigned int s = 0 ; s < current_cell.sphereId.size() ; ++s)
- {
- id = current_cell.sphereId[s];
- if (id != sphereId && sphere[id].R > 0.0)
- {
- if ((distance_centre_spheres(S,sphere[id]) - (S.R + sphere[id].R)) < -max_overlap_rate * rmin) { failure = 128; break; }
- }
- }
- }
- }
- }
-
- }
+ failure = check_overlaps(S, sphereId);
+ else continue;
if (!failure)
{
@@ -707,6 +845,33 @@
return 0;
}
+unsigned int SpherePadder::check_overlaps(Sphere & S, unsigned int excludedId)
+{
+ unsigned int id;
+ Cell current_cell;
+
+ partition.locateCellOf(S.x,S.y,S.z);
+ for (unsigned int i = partition.i_down() ; i <= partition.i_up() ; ++i)
+ {
+ for (unsigned int j = partition.j_down() ; j <= partition.j_up() ; ++j)
+ {
+ for (unsigned int k = partition.k_down() ; k <= partition.k_up() ; ++k)
+ {
+ current_cell = partition.get_cell(i,j,k);
+ for (unsigned int s = 0 ; s < current_cell.sphereId.size() ; ++s)
+ {
+ id = current_cell.sphereId[s];
+ if (id != excludedId && sphere[id].R > 0.0)
+ {
+ if (distance_spheres(S,sphere[id]) < (-max_overlap_rate * rmin)) { return 128; }
+ }
+ }
+ }
+ }
+ }
+ return 0;
+}
+
double SpherePadder::distance_vector3 (double V1[],double V2[])
{
double D[3];
@@ -723,15 +888,18 @@
C1[0] = sphere[s1].x ; C1[1] = sphere[s1].y ; C1[2] = sphere[s1].z ; R1 = sphere[s1].R;
C2[0] = sphere[s2].x ; C2[1] = sphere[s2].y ; C2[2] = sphere[s2].z ; R2 = sphere[s2].R;
C3[0] = sphere[s3].x ; C3[1] = sphere[s3].y ; C3[2] = sphere[s3].z ; R3 = sphere[s3].R;
- C4[0] = sphere[s4].x ; C4[1] = sphere[s4].y ; C4[2] = sphere[s4].z ; R4 = sphere[s4].R;
+ C4[0] = sphere[s4].x ; C4[1] = sphere[s4].y ; C4[2] = sphere[s4].z ; R4 = sphere[s4].R;
+
+ // rougth estimate of the position
+ // TODO
- unsigned int fail_det = 1;
- unsigned int fail_delta = 2;
- unsigned int fail_radius = 4;
- unsigned int fail_overlap = 8;
- unsigned int fail_gap = 16;
- unsigned int fail_radius_range = 32;
- unsigned int fail_NaN = 64;
+ unsigned int fail_det = 0x01;
+ unsigned int fail_delta = 0x02;
+ unsigned int fail_radius = 0x04;
+ unsigned int fail_overlap = 0x08;
+ unsigned int fail_gap = 0x10;
+ unsigned int fail_radius_range = 0x20;
+ unsigned int fail_NaN = 0x40;
// (x-x1)^2 + (y-y1)^2 + (z-z1)^2 = (r+r1)^2 (1)
// (x-x2)^2 + (y-y2)^2 + (z-z2)^2 = (r+r2)^2 (2)
@@ -760,9 +928,9 @@
double eee = (C1[0]*C1[0] + C1[1]*C1[1] + C1[2]*C1[2] - R1*R1) - (C4[0]*C4[0] + C4[1]*C4[1] + C4[2]*C4[2] - R4*R4);
// compute the determinant of matrix A (system AX = B)
- // [a , b, c]; [x] [e - d*r]
- // A = [aa ,bb ,cc ]; X=[y] B=[ee - dd*r]
- // [aaa,bbb,ccc]; [z] [eee-ddd*r]
+ // [a , b, c]; [x] [e - d*r]
+ // A = [aa ,bb ,cc ]; X=[y] B(r)=[ee - dd*r]
+ // [aaa,bbb,ccc]; [z] [eee-ddd*r]
double DET = a*(bb*ccc-bbb*cc) - aa*(b*ccc-bbb*c) + aaa*(b*cc-bb*c);
double R = 0.0;
@@ -784,6 +952,7 @@
double a32 = -(a*bbb-aaa*b) * inv_DET;
double a33 = (a*bb-aa*b) * inv_DET;
+ // A^-1 B(r)
double xa = -(a11*d + a12*dd + a13*ddd);
double xb = (a11*e + a12*ee + a13*eee);
@@ -792,7 +961,7 @@
double za = -(a31*d + a32*dd + a33*ddd);
double zb = (a31*e + a32*ee + a33*eee);
-
+
// Replace x, y and z in Equation (1) and solve the second order equation A*r^2 + B*r + C = 0
double A = xa*xa + ya*ya + za*za - 1.0;
@@ -802,9 +971,9 @@
double DELTA = B*B - 4.0*A*C;
double RR1,RR2;
- if (DELTA >= 0.0)
+ if (DELTA >= 0.0) // && A != 0.0
{
- double inv_2A = 1.0 / (2.0*A);
+ double inv_2A = 0.5/A;
double sqrt_DELTA = sqrt(DELTA);
RR1 = (-B + sqrt_DELTA) * inv_2A;
RR2 = (-B - sqrt_DELTA) * inv_2A;
@@ -822,6 +991,11 @@
centre[2] = za * R + zb;
}
else return fail_det;
+
+ S.x = centre[0];
+ S.y = centre[1];
+ S.z = centre[2];
+ S.R = R;
// FIXME do not use sqrt to speed up... (squared_distance_vector3)
// for the moment it is not critical
@@ -832,40 +1006,42 @@
double distance3 = distance_vector3 (centre,C3) - (R + R3);
double distance4 = distance_vector3 (centre,C4) - (R + R4);
- double half_distance_rate = -0.5 * max_overlap_rate;
- if ( ( distance1 < half_distance_rate * (R + R1))
- || ( distance2 < half_distance_rate * (R + R2))
- || ( distance3 < half_distance_rate * (R + R3))
- || ( distance4 < half_distance_rate * (R + R4)) )
+// double half_distance_rate = -0.5 * max_overlap_rate;
+// if ( ( distance1 < half_distance_rate * (R + R1))
+// || ( distance2 < half_distance_rate * (R + R2))
+// || ( distance3 < half_distance_rate * (R + R3))
+// || ( distance4 < half_distance_rate * (R + R4)) )
+// { return fail_overlap; }
+
+ double max_overlap = -max_overlap_rate * rmin;
+ if ( ( distance1 < max_overlap)
+ || ( distance2 < max_overlap)
+ || ( distance3 < max_overlap)
+ || ( distance4 < max_overlap) )
{ return fail_overlap; }
// The gap between spheres must not be too large
double distance_max = max_overlap_rate * rmin;
- if ( ( distance1 > distance_max)
- || ( distance2 > distance_max)
- || ( distance3 > distance_max)
- || ( distance4 > distance_max) )
- { return fail_gap; }
+ if ( ( distance1 > distance_max)
+ || ( distance2 > distance_max)
+ || ( distance3 > distance_max)
+ || ( distance4 > distance_max) )
+ { return fail_gap; }
-
// Check if there is NaN
- if ( (centre[0] != centre[0])
- || (centre[1] != centre[1])
- || (centre[2] != centre[2])
- || (R != R))
+ if ( (centre[0] != centre[0])
+ || (centre[1] != centre[1])
+ || (centre[2] != centre[2])
+ || (R != R))
{ return fail_NaN; }
-
- S.x = centre[0];
- S.y = centre[1];
- S.z = centre[2];
- S.R = R;
-
+
+
+ //S.R = R;
return 0;
}
-
- void SpherePadder::place_virtual_spheres()
+ void SpherePadder::place_virtual_spheres() // WARNING there is a bug in this function (TODO include modifications of Philippe Marin)
{
BEGIN_FUNCTION("Place virtual spheres...");
@@ -958,6 +1134,7 @@
if ( k < 0) k = -k;
if ( scalar_product < 0) k = -k;
// TODO : il suffit de tester la valeur de mesh->face[f].normal_swap
+ // mais pour le moment cette routine n'est pas critique (suffisemment rapide)
// First virtual sphere
@@ -1101,8 +1278,8 @@
S.R = R;
S.type = INSERTED_BY_USER;
+ unsigned int inseredId = sphere.size();
sphere.push_back(S);
- unsigned int inseredId = sphere.size()-1;
partition.locateCellOf(x-R, y-R, z-R);
unsigned int iCellMin = partition.i_down();
@@ -1118,19 +1295,20 @@
unsigned int id;
Cell current_cell;
- for (unsigned int i = iCellMin ; i <= iCellMax ; ++i)
+ for (unsigned int i = iCellMin ; i <= iCellMax ; ++i /*unsigned int i = 0 ; i <= partition.isize-1 ; ++i*/)
{
- for (unsigned int j = jCellMin ; j <= jCellMax ; ++j)
+ for (unsigned int j = jCellMin ; j <= jCellMax ; ++j /*unsigned int j = 0 ; j <= partition.jsize-1 ; ++j*/)
{
- for (unsigned int k = kCellMin ; k <= kCellMax ; ++k)
+ for (unsigned int k = kCellMin ; k <= kCellMax ; ++k /*unsigned int k = 0 ; k <= partition.ksize-1 ; ++k*/)
{
current_cell = partition.get_cell(i,j,k);
+ partition.add_in_cell(inseredId,i,j,k); // FIXME could be set only in iCellMin+1 to iCellMax-1 ...
for (unsigned int s = 0 ; s < current_cell.sphereId.size() ; ++s)
{
id = current_cell.sphereId[s];
- if (sphere[id].type != INSERTED_BY_USER && sphere[id].R > 0.0)
+ if (sphere[id].type != INSERTED_BY_USER && sphere[id].type != VIRTUAL && sphere[id].R > 0.0)
{
- if (distance_spheres(inseredId,id) < distance_max )
+ if (distance_spheres(inseredId,id) < distance_max)
{
sphere[id].R = 0.0;
}
@@ -1139,7 +1317,29 @@
}
}
}
-
+
+
+// for (/*unsigned int i = iCellMin ; i <= iCellMax ; ++i*/unsigned int i = 0 ; i <= partition.isize-1 ; ++i)
+// {
+// for (/*unsigned int j = jCellMin ; j <= jCellMax ; ++j*/unsigned int j = 0 ; j <= partition.jsize-1 ; ++j)
+// {
+// for (/*unsigned int k = kCellMin ; k <= kCellMax ; ++k*/unsigned int k = 0 ; k <= partition.ksize-1 ; ++k)
+// {
+//
+// current_cell = partition.get_cell(i,j,k);
+// for (unsigned int s = 0 ; s < current_cell.sphereId.size() ; ++s)
+// {
+//
+// if (current_cell.sphereId[s] == inseredId)
+// cout << "Cell " << i << ", " << j << ", " << k << " OK -> " << current_cell.sphereId[s] << endl;
+//
+// }
+// }
+// }
+// }
+
+
+
}
@@ -1170,12 +1370,21 @@
id = current_cell.sphereId[s];
if (id != sphereId && sphere[id].R > 0.0)
{
+ // Virtual spheres and spheres added by user cannot be modified here
+// unsigned int opt = 0x00;
+// if (sphere[sphereId].type == VIRTUAL || sphere[sphereId].type == INSERTED_BY_USER) opt |= 0x01;
+// if (sphere[id].type == VIRTUAL || sphere[id].type == INSERTED_BY_USER) opt |= 0x02;
+// if (opt == 0x03) continue;
+ double inv_sumR = 1.0 / (sphere[sphereId].R + sphere[id].R);
while ( (distance = distance_spheres(sphereId,id)) < distance_max )
{
- kfact = 1.0 + distance / (sphere[sphereId].R + sphere[id].R);
- sphere[sphereId].R *= kfact;
- sphere[id].R *= kfact;
+ kfact = 1.0 + distance * inv_sumR;
+ //if (opt | 0x01) // test a revoir
+ sphere[sphereId].R *= kfact;
+ //if (opt | 0x02)
+ sphere[id].R *= kfact;
}
+
if (sphere[id].R < rmin) sphere[id].R = 0.0;
if (sphere[sphereId].R < rmin) sphere[sphereId].R = 0.0;
}
@@ -1190,11 +1399,9 @@
}
-// FIXME attention aux interpenetrations avec rayons nuls
- void SpherePadder::detect_overlap ()
+void SpherePadder::detect_overlap ()
{
BEGIN_FUNCTION("Detect overlap");
- cerr << endl;
Sphere S1,S2;
unsigned int nb_overlap = 0;
double d;
@@ -1202,7 +1409,7 @@
for (unsigned int n = 0 ; n < sphere.size()-1 ; ++n)
{
S1 = sphere[ n ];
- if ( S1.type != VIRTUAL && S1.R > rmin)
+ if ( S1.type != VIRTUAL && S1.R > rmin )
{
for (unsigned int t = n+1 ; t < sphere.size() ; ++t)
{
@@ -1215,17 +1422,18 @@
{
cerr << "Overlap!" << endl;
partition.locateCellOf(S1.x,S1.y,S1.z);
- cerr << "ijk (S1) = " << partition.current_i << ", " << partition.current_j << ", " << partition.current_k << endl;
+ cerr << " Spheres = " << n << ", " << t << endl;
+ cerr << " Cell ijk (S1) = " << partition.current_i << ", " << partition.current_j << ", " << partition.current_k << endl;
partition.locateCellOf(S2.x,S2.y,S2.z);
- cerr << "ijk (S2) = " << partition.current_i << ", " << partition.current_j << ", " << partition.current_k << endl;
+ cerr << " Cell ijk (S2) = " << partition.current_i << ", " << partition.current_j << ", " << partition.current_k << endl;
cerr << " Types = " << S1.type << ", " << S2.type << endl;
cerr << " Radii = " << S1.R << ", " << S2.R << endl;
cerr << " pos S1 = " << S1.x << ", " << S1.y << ", " << S1.z << endl;
cerr << " pos S2 = " << S2.x << ", " << S2.y << ", " << S2.z << endl;
cerr << " Distance / rmin = " << d/rmin << endl;
- cerr << " rmin = " << rmin << endl;
- sphere[n].type = sphere[t].type = AT_TETRA_VERTEX;
- sphere[n].R = sphere[t].R = 0.0;
+ //cerr << " rmin = " << rmin << endl;
+ //sphere[n].type = sphere[t].type = AT_TETRA_VERTEX;
+ //sphere[n].R = sphere[t].R = 0.0;
++(nb_overlap);
}
}
@@ -1275,7 +1483,6 @@
volume_portion4 = spherical_triangle(point4, point1, point2, point3);
return (volume_portion1 + volume_portion2 + volume_portion3 + volume_portion4);
-
}
// point : x,y,z,R
@@ -1319,9 +1526,72 @@
}
+// ============================
+// ANALISYS
+// ============================
+void SpherePadder::check_inProbe(unsigned int i)
+{
+ if (!probeIsDefined) return;
+ if (sphere[i].R <= 0.0) return;
+
+ double dx = sphere[i].x - xProbe;
+ double dy = sphere[i].y - yProbe;
+ double dz = sphere[i].z - zProbe;
+ double squared_distance = dx*dx + dy*dy + dz*dz;
+ if (squared_distance + sphere[i].R*sphere[i].R < RProbe * RProbe) sphereInProbe.push_back(i);
+}
+void SpherePadder::add_spherical_probe(double Rfact)
+{
+ double xmin,xmax;
+ double ymin,ymax;
+ double zmin,zmax;
+ double R;
+
+ xmin=xmax=mesh->node[0].x;
+ ymin=ymax=mesh->node[0].y;
+ zmin=zmax=mesh->node[0].z;
+ for (unsigned int i = 1 ; i < mesh->node.size() ; ++i)
+ {
+ xmin = (xmin > mesh->node[i].x) ? mesh->node[i].x : xmin;
+ xmax = (xmax < mesh->node[i].x) ? mesh->node[i].x : xmax;
+ ymin = (ymin > mesh->node[i].y) ? mesh->node[i].y : ymin;
+ ymax = (ymax < mesh->node[i].y) ? mesh->node[i].y : ymax;
+ zmin = (zmin > mesh->node[i].z) ? mesh->node[i].z : zmin;
+ zmax = (zmax < mesh->node[i].z) ? mesh->node[i].z : zmax;
+ }
+ xProbe = 0.5 * (xmin + xmax);
+ yProbe = 0.5 * (ymin + ymax);
+ zProbe = 0.5 * (zmin + zmax);
+ RProbe = (xmax - xmin);
+ RProbe = (RProbe < (R = ymax - ymin)) ? RProbe : R;
+ RProbe = (RProbe < (R = zmax - zmin)) ? RProbe : R;
+ RProbe *= 0.5 * Rfact;
+ probeIsDefined = true;
+ compacity_file.open("compacity.dat");
+}
+double SpherePadder::compacity_in_probe(unsigned int ninsered)
+{
+ if (!probeIsDefined) return -1.0;
+
+ double dr = 0.0;
+ double Vs = 0.0;
+ double fact = 1.333333333 * M_PI;
+ for (unsigned int i = 1 ; i < sphereInProbe.size() ; ++i)
+ {
+ Vs += fact * sphere[i].R * sphere[i].R * sphere[i].R;
+ }
+
+ // TODO ameliorer cette fonction en prennant en compte les spheres coupees...
+
+ dr = Vs / (fact*RProbe*RProbe*RProbe);
+ compacity_file << ninsered << ' ' << dr << endl;
+ return (dr);
+}
+
+
Modified: trunk/extra/SpherePadder/SpherePadder.hpp
===================================================================
--- trunk/extra/SpherePadder/SpherePadder.hpp 2009-04-06 10:26:50 UTC (rev 1746)
+++ trunk/extra/SpherePadder/SpherePadder.hpp 2009-04-06 15:42:25 UTC (rev 1747)
@@ -16,7 +16,7 @@
#include <time.h>
#include <set>
-# define BEGIN_FUNCTION(arg) if (trace_functions) cerr << (arg) << "... " << flush
+# define BEGIN_FUNCTION(arg) if (trace_functions) cerr << (arg) << "..." << endl << flush
# define END_FUNCTION if (trace_functions) cerr << "Done\n" << flush
enum SphereType {AT_NODE, AT_SEGMENT, AT_FACE, AT_TETRA_CENTER, AT_TETRA_VERTEX, INSERTED_BY_USER, FROM_TRIANGULATION, VIRTUAL};
@@ -25,7 +25,7 @@
{
double x,y,z,R;
SphereType type;
- unsigned int tetraOwner;
+ unsigned int tetraOwner; // FIXME can be removed ??
};
struct Neighbor
@@ -62,12 +62,16 @@
protected:
vector<vector<unsigned int> > combination;// FIXME long ?
- SpherePackingTriangulation triangulation;// TODO use Delaunay Triangulation to avoid flat tetrahedra
+ SpherePackingTriangulation triangulation;
vector<tetra_porosity> tetra_porosities;
double distance_spheres (unsigned int i, unsigned int j);// FIXME long ?
+ double distance_spheres (Sphere& S1,Sphere& S2);
+ double squared_distance_centre_spheres(unsigned int i, unsigned int j);
double distance_centre_spheres(Sphere& S1, Sphere& S2);
double distance_vector3 (double V1[],double V2[]);
+ void build_sorted_list_of_neighbors(unsigned int sphereId, vector<neighbor_with_distance> & neighbor);
+ void build_sorted_list_of_neighbors(Sphere & S, vector<neighbor_with_distance> & neighbor);
double spherical_triangle (double point1[],double point2[],double point3[],double point4[]);
double solid_volume_of_tetrahedron(Sphere& S1, Sphere& S2, Sphere& S3, Sphere& S4);
void place_at_nodes ();
@@ -76,20 +80,20 @@
void place_at_tetra_centers ();
void place_at_tetra_vertexes ();
void cancel_overlaps ();
+ unsigned int iter_densify(unsigned int nb_check = 100);
+ void repack_null_radii();
-
- //
+ // some key functions
unsigned int place_fifth_sphere(unsigned int s1, unsigned int s2, unsigned int s3, unsigned int s4, Sphere& S);// FIXME long ?
- unsigned int place_sphere_4contacts (unsigned int sphereId, unsigned int nb_combi_max = 30);// FIXME long ?
-
- // Check functions
- void detect_overlap ();
-
+ unsigned int place_sphere_4contacts (unsigned int sphereId, unsigned int nb_combi_max = 15);// FIXME long ?
+ unsigned int check_overlaps(Sphere & S, unsigned int excludedId);
+
double rmin,rmax,rmoy,dr;
- double ratio;
+ double ratio; // rmax/rmin
double max_overlap_rate;
unsigned int n1,n2,n3,n4,n5,n_densify;
unsigned int nb_iter_max;
+ bool RadiusDataIsOK,RadiusIsSet;
TetraMesh * mesh;
vector <Sphere> sphere;
@@ -111,17 +115,28 @@
public:
bool meshIsPlugged;
-
+
+ void setRadiusRatio(double r);
+ void setRadiusRange(double min, double max);
+ void setMaxOverlapRate(double r) { max_overlap_rate = fabs(r); }
+ void setVirtualRadiusFactor(double f) {virtual_radius_factor = fabs(f);}
+ void setTargetSolidFraction(double sf)
+ {
+ // TODO
+ }
+
void plugTetraMesh (TetraMesh * mesh);
void save_mgpost (const char* name);
- void save_Rxyz (const char* name);
+ void save_tri_mgpost (const char* name);
+ void save_Rxyz (const char* name);
double virtual_radius_factor;
SpherePadder();
-
- // Pading functions
+ // Check functions (to debug)
+ void detect_overlap ();
+
//! \brief 5-step padding (for details see Jerier et al.)
void pad_5 ();
Modified: trunk/extra/SpherePadder/main.cpp
===================================================================
--- trunk/extra/SpherePadder/main.cpp 2009-04-06 10:26:50 UTC (rev 1746)
+++ trunk/extra/SpherePadder/main.cpp 2009-04-06 15:42:25 UTC (rev 1747)
@@ -10,8 +10,6 @@
#include "SpherePadder.hpp"
-#define DEBUG
-
unsigned int mesh_format;
vector <unsigned int> output_format;
char mesh_file_name[100];
@@ -20,7 +18,7 @@
int main()
{
-#ifdef DEBUG
+#if 1
TetraMesh * mesh = new TetraMesh();
//mesh->read_gmsh("meshes/cube1194.msh");
@@ -28,14 +26,19 @@
//mesh->write_surface_MGP ("cube.mgp");
SpherePadder * padder = new SpherePadder();
+
padder->plugTetraMesh(mesh);
- //padder->add_spherical_probe(0.7);
-
+ padder->setRadiusRatio(4.0);
+ padder->setMaxOverlapRate(1.0e-4);
+ padder->setVirtualRadiusFactor(100.0);
+
padder->pad_5();
- padder->insert_sphere(0.5,0.5,0.5,0.2);
padder->place_virtual_spheres();
+ padder->insert_sphere(0.5,0.5,0.5,0.4);
padder->densify();
-
+ padder->detect_overlap ();
+
+ //padder->save_tri_mgpost("triangulation.mgp");
padder->save_mgpost("mgp.out.001");
padder->save_Rxyz("spheres.Rxyz");
return 0;
