------------------------------------------------------------
revno: 2991
committer: Bruno Chareyre <bruno.chare...@hmg.inpg.fr>
branch nick: trunk
timestamp: Fri 2012-01-13 12:24:44 +0100
message:
- make deformation analysis possible for arbitrary simulations (no longer
depends on TriaxialEngine for defining the bounds)
modified:
lib/triangulation/KinematicLocalisationAnalyser.cpp
lib/triangulation/KinematicLocalisationAnalyser.hpp
pkg/dem/MicroMacroAnalyser.cpp
pkg/dem/MicroMacroAnalyser.hpp
pkg/dem/TesselationWrapper.cpp
pkg/dem/TesselationWrapper.hpp
--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to
https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription
=== modified file 'lib/triangulation/KinematicLocalisationAnalyser.cpp'
--- lib/triangulation/KinematicLocalisationAnalyser.cpp 2012-01-11 21:30:18 +0000
+++ lib/triangulation/KinematicLocalisationAnalyser.cpp 2012-01-13 11:24:44 +0000
@@ -23,7 +23,6 @@
{
int n_debug = 0;
-//Usefull fonction to convert int to string (define it elsewhere)
std::string _itoa(int i)
{
std::ostringstream buffer;
@@ -95,23 +94,17 @@
file_number_1 = n1;
file_number_0 = n0;
base_file_name = string(base_name);
-
consecutive = ((n1-n0)==1);
bz2 = usebz2;
sphere_discretisation = SPHERE_DISCRETISATION;
linear_discretisation = LINEAR_DISCRETISATION;
TS1 = new(TriaxialState);
TS0 = new(TriaxialState);
-
-// char buffer [50];
std::ostringstream file_name1, file_name0;
file_name1 << (string)(base_file_name) << n1;
file_name0 << (string)(base_file_name) << n0;
- //cout << "file names : " << file_name0.str().c_str() << ", " << file_name1.str().c_str() << endl;
-
TS1->from_file(file_name1.str().c_str(), bz2);
TS0->from_file(file_name0.str().c_str(), bz2);
-
Delta_epsilon(3,3) = TS1->eps3 - TS0->eps3;
Delta_epsilon(1,1) = TS1->eps1 - TS0->eps1;
Delta_epsilon(2,2) = TS1->eps2 - TS0->eps2;
@@ -122,14 +115,10 @@
base_file_name = name;
}
-
bool KinematicLocalisationAnalyser::SetFileNumbers(int n0, int n1)
{
- //TriaxialState* TS3;
- //string file_name;
bool bf0 = false;
bool bf1 = false;
-// char buffer [50];
if (file_number_0 != n0) {
if (file_number_1 != n0) {
//file_name = base_file_name + n0;
@@ -155,7 +144,6 @@
return (bf0 && bf1);
}
-
void KinematicLocalisationAnalyser::SetConsecutive(bool t)
{
consecutive = t;
@@ -167,7 +155,6 @@
TS1->NO_ZERO_ID = t;
}
-
void KinematicLocalisationAnalyser::SwitchStates(void)
{
TriaxialState* TStemp = TS0;
@@ -213,22 +200,22 @@
vtk.begin_vertices();
RTriangulation::Finite_vertices_iterator V_it = Tri.finite_vertices_begin();
+ bool beginWithFictious = V_it->info().isFictious;
for (; V_it != Tri.finite_vertices_end(); ++V_it) if (!V_it->info().isFictious) vtk.file << V_it->point().point() << endl;
vtk.end_vertices();
vtk.begin_cells();
Finite_cells_iterator cell = Tri.finite_cells_begin();
- //FIXME : Preconditions : scene has 6 boxes with ids 0-5, and bodies are listed in ascending ids order
+ //FIXME : Preconditions : the fictious bounds are first in the list
for (; cell != Tri.finite_cells_end(); ++cell) {
if (!cell->info().isFictious) vtk.write_cell(
- cell->vertex(0)->info().id()-6,
- cell->vertex(1)->info().id()-6,
- cell->vertex(2)->info().id()-6,
- cell->vertex(3)->info().id()-6
+ cell->vertex(0)->info().id()- (beginWithFictious?n_fictious_vertices:0),
+ cell->vertex(1)->info().id()- (beginWithFictious?n_fictious_vertices:0),
+ cell->vertex(2)->info().id()- (beginWithFictious?n_fictious_vertices:0),
+ cell->vertex(3)->info().id()- (beginWithFictious?n_fictious_vertices:0)
);
}
vtk.end_cells();
-
vtk.begin_data("Strain_matrix", POINT_DATA, TENSORS, FLOAT);
V_it = Tri.finite_vertices_begin();
for (; V_it != Tri.finite_vertices_end(); ++V_it) {
@@ -237,7 +224,6 @@
vtk.file << ParticleDeformation[V_it->info().id()] << endl;}
}
vtk.end_data();
-
vtk.begin_data("Strain_deviator", POINT_DATA, SCALARS, FLOAT);
V_it = Tri.finite_vertices_begin();
for (; V_it != Tri.finite_vertices_end(); ++V_it) {
@@ -248,44 +234,7 @@
//vtk.write_data((float) epsilon.Deviatoric()(1,1)-epsilon.Deviatoric()(0,0));}
}
vtk.end_data();
-
return true;
-
- /*
- ofstream output_file(output_file_name);
- if (!output_file.is_open()) {
- cerr << "Error opening files";
- return false;
- }
- ComputeParticlesDeformation();
-
- Tesselation& Tes = TS1->tesselation();
- RTriangulation& Tri = Tes.Triangulation();
-
- output_file << Tri.number_of_vertices()<<endl;
- for (RTriangulation::Finite_vertices_iterator V_it =
- Tri.finite_vertices_begin(); V_it != Tri.finite_vertices_end(); V_it++)
- output_file<<V_it->info().id()<<" "<<V_it->point()<<endl;
-
- output_file << Tri.number_of_finite_cells()<<endl;
- Finite_cells_iterator cell = Tri.finite_cells_begin();
- Finite_cells_iterator cell0 = Tri.finite_cells_end();
- for (; cell != cell0; cell++) {
- for (unsigned int index=0; index<4; index++) output_file << cell->vertex(index)->info().id()<<" " ;
- output_file<<endl;
- }
-
- for (RTriangulation::Finite_vertices_iterator V_it =
- Tri.finite_vertices_begin(); V_it != Tri.finite_vertices_end(); V_it++) {
- Tenseur_sym3 epsilon(ParticleDeformation[V_it->info().id()]); // partie symetrique
- double dev = (double) epsilon.Deviatoric().Norme();
- output_file<<V_it->info().id()<<endl<<ParticleDeformation[V_it->info().id()]<<dev<<endl;
-
- }
- */
- return 1;
-
-
}
bool KinematicLocalisationAnalyser::DistribsToFile(const char* output_file_name)
@@ -293,8 +242,7 @@
ofstream output_file(output_file_name);
if (!output_file.is_open()) {
cerr << "Error opening files";
- return false;
- }
+ return false;}
output_file << "sym_grad_u_total_g (wrong averaged strain):"<< endl << Tenseur_sym3(grad_u_total_g) << endl;
output_file << "Total volume = " << v_total << ", grad_u = " << endl << grad_u_total << endl << "sym_grad_u (true average strain): " << endl << Tenseur_sym3(grad_u_total) << endl;
@@ -354,8 +302,6 @@
}
}
NormalDisplacementDistributionToFile(edges, output_file);
-
-
output_file.close();
return true;
}
@@ -368,8 +314,7 @@
if (state.inside((*cit)->grain1->sphere.point()) && state.inside((*cit)->grain2->sphere.point()))
nc1 += 2;
else if (state.inside((*cit)->grain1->sphere.point()) || state.inside((*cit)->grain2->sphere.point()))
- ++nc1;
- }
+ ++nc1;}
return nc1;
}
@@ -465,8 +410,7 @@
state.inside((*cit)->grain2->sphere.point())) {
v = (*cit)->normal;
for (int i=1; i<4; i++) for (int j=3; j>=i; j--)
- Tens(i,j) += v[i-1]*v[j-1];
- }
+ Tens(i,j) += v[i-1]*v[j-1];}
}
Tens /= Filtered_contacts(state);
return Tens;
@@ -566,7 +510,6 @@
}
output_file << endl;
return output_file;
-
}
ofstream& KinematicLocalisationAnalyser::
@@ -624,13 +567,11 @@
output_file << row[i].first << " " << row[i].second << endl;
cerr << row[i].first << " " << row[i].second << endl;
}
-
output_file << endl;
return output_file;
}
-
ofstream& KinematicLocalisationAnalyser::
AllNeighborDistributionToFile(ofstream& output_file)
{
@@ -689,7 +630,6 @@
output_file << row[i].first << " " << row[i].second << endl;
cerr << row[i].first << " " << row[i].second << endl;
}
-
output_file << endl;
return output_file;
}
@@ -703,12 +643,6 @@
// Real DZ = 1.0/sphere_discretisation;
long Nc0 = TS0->contacts.size();
long Nc1 = TS1->contacts.size();
-
-// long nv1=0;
-// long nv2=0;
-// long nv3=0;
-// long ng1=0;
-// long ng2=0;
n_persistent = 0; n_new = 0; n_lost = 0;
long lost_in_state0 = 0;
@@ -766,57 +700,8 @@
++n_new;
}
}
- /*
- RGrid1D table;
-
- for ( Edge_iterator ed_it = T.edges_begin(); ed_it != T.edges_end(); ed_it++ )
- {
- if ( !T.is_infinite ( *ed_it ) )
- {
- s = T.segment ( *ed_it );
- if ( ( *TS1 ).inside ( s.source() ) && ( *TS1 ).inside ( s.target() ) )
- {
- v = s.to_vector();
- row[ ( int ) ( abs ( v.z() /sqrt ( s.squared_length() ) ) /DZ ) ].second += 2;
- nv1 += 2;
- }
- else
- {
- if ( ( *TS1 ).inside ( s.source() ) || ( *TS1 ).inside ( s.target() ) )
- {
- v = s.to_vector();
- row[ ( int ) ( abs ( v.z() /sqrt ( s.squared_length() ) ) /DZ ) ].second += 1;
- ++nv1;
- }
- else ++nv2;
- }
- }
- else ++nv3;
- }
-
- Real normalize = 1.0/ ( ng1*4*DZ*3.141592653 );
- for ( int i = 0; i < sphere_discretisation; ++i ) row[i].second *= normalize;
-
- output_file << "#Neighbors distribution" << endl << "(filter dist. = " << ( *TS1 ).filter_distance
- << ", "<< nv1 << " neighbors + " << nv2 << " excluded + "
- << nv3 << " infinite, for "<< ng1 <<"/"<< ( ng1+ng2 ) << " grains)" << endl;
- output_file << "max_nz number_of_neighbors" << endl;
- cerr << "#Neighbors distribution" << endl << "(filter dist. = " << ( *TS1 ).filter_distance
- << ", "<< nv1 << " neighbors + " << nv2 << " excluded + "
- << nv3 << " infinite, for "<< ng1 <<"/"<< ( ng1+ng2 ) << " grains)" << endl;
- cerr << "mean_nz number_of_neighbors" << endl;
- for ( int i = 0; i < sphere_discretisation; ++i )
- {
- output_file << row[i].first << " " << row[i].second << endl;
- cerr << row[i].first << " " << row[i].second << endl;
- }
-
- output_file << endl;
- return output_file;*/
}
-
-
void KinematicLocalisationAnalyser::SetDisplacementIncrements(void)
{
TriaxialState::GrainIterator gend = TS1->grains_end();
@@ -835,7 +720,6 @@
{
Vecteur v(0.f, 0.f, 0.f);
int id;// ident. de la particule
-
Vecteur fixedPoint = 0.5*((TS0->box.base-CGAL::ORIGIN)+(TS0->box.sommet-CGAL::ORIGIN));
for (int i=0; i<4; i++) {
// char msg [256];
@@ -848,24 +732,14 @@
meanFieldDisp[1]*Delta_epsilon(1,1),
meanFieldDisp[2]*Delta_epsilon(2,2));
} else meanFieldDisp=Vecteur(0,0,0);
- if (consecutive)
- v = v + TS1->grain(id).translation-meanFieldDisp;
- else {
- v = v + (TS1->grain(id).sphere.point() - TS0->grain(id).sphere.point()-meanFieldDisp);
-
- }
-
- //for tests with affine displacement field
- //if ((TS1->grain(id).sphere.point().y()+TS1->grain(id).sphere.point().z())>0.035)//a discontinuity
- //v = v + Vecteur(0, 0.01*TS1->grain(id).sphere.point().x(), 0);
+ if (consecutive) v = v + TS1->grain(id).translation-meanFieldDisp;
+ else v = v + (TS1->grain(id).sphere.point() - TS0->grain(id).sphere.point()-meanFieldDisp);
}
}
- v = v*0.333333333333333333333333;
+ v = v*0.333333333333;
return v;
}
-
-
void KinematicLocalisationAnalyser::Grad_u(Finite_cells_iterator cell, int facet, Vecteur &V, Tenseur3& T)
{
Vecteur S = cross_product((cell->vertex(l_vertices[facet][1])->point())
@@ -875,10 +749,9 @@
Somme(T, V, S);
}
-
void KinematicLocalisationAnalyser::Grad_u(Finite_cells_iterator cell,
Tenseur3& T, bool vol_divide)// Calcule le gradient de d�p.
-{
+{
T.reset();
Vecteur v;
for (int facet=0; facet<4; facet++) {
@@ -910,12 +783,12 @@
ParticleDeformation.resize(Tes.Max_id() + 1);
}
//reset volumes and tensors of each particle
- n_real_vertices = 0;
+ n_real_vertices = 0; n_fictious_vertices=0;
for (RTriangulation::Finite_vertices_iterator V_it=Tri.finite_vertices_begin(); V_it != Tri.finite_vertices_end(); V_it++) {
//cerr << V_it->info().id() << endl;
V_it->info().v() =0;//WARNING : this will erase previous values if some have been computed
ParticleDeformation[V_it->info().id()]=NULL_TENSEUR3;
- if (!V_it->info().isFictious) ++n_real_vertices;
+ if (!V_it->info().isFictious) ++n_real_vertices; else ++n_fictious_vertices;
}
Finite_cells_iterator cell = Tri.finite_cells_begin();
Finite_cells_iterator cell0 = Tri.finite_cells_end();
=== modified file 'lib/triangulation/KinematicLocalisationAnalyser.hpp'
--- lib/triangulation/KinematicLocalisationAnalyser.hpp 2012-01-11 21:30:18 +0000
+++ lib/triangulation/KinematicLocalisationAnalyser.hpp 2012-01-13 11:24:44 +0000
@@ -113,7 +113,7 @@
Real v_solid_total;//solid volume in the box
Real v_total;//volume of the box
Real v_total_g;//summed volumes of extended grain cells
- long n_persistent, n_new, n_lost, n_real_cells, n_real_vertices;
+ long n_persistent, n_new, n_lost, n_real_cells, n_real_vertices, n_fictious_vertices;
=== modified file 'pkg/dem/MicroMacroAnalyser.cpp'
--- pkg/dem/MicroMacroAnalyser.cpp 2010-11-19 12:30:08 +0000
+++ pkg/dem/MicroMacroAnalyser.cpp 2012-01-13 11:24:44 +0000
@@ -113,11 +113,14 @@
BodyContainer::iterator bi = biBegin;
Ng = 0;
vector<Body::id_t> fictiousVtx;
+ bool fictiousFirst=false;
for (; bi!=biEnd ; ++bi) {
const Body::id_t Idg = (*bi)->getId();
TS.grains[Idg].id = Idg;
if (!(*bi)->isDynamic()) {
- TS.grains[Idg].isSphere = false; fictiousVtx.push_back(Idg);}
+ if (!nonSphereAsFictious) continue;
+ TS.grains[Idg].isSphere = false; fictiousVtx.push_back(Idg);
+ if (bi==biBegin) fictiousFirst=true;}
else {//then it is a sphere (not a wall)
++Ng;
const Sphere* s = YADE_CAST<Sphere*> ((*bi)->shape.get());
@@ -141,9 +144,11 @@
}
TS.mean_radius /= Ng;//rayon moyen
LOG_INFO(" loaded : " << Ng << " grains with mean radius = " << TS.mean_radius);
- if (fictiousVtx.size()>=6){//boxes found, simulate them with big spheres
+ Real FAR = 1e4;
+ if (fictiousVtx.size()==0) {
+ TS.grains.resize(Ng+6); for (int fv=Ng;fv<Ng+6;fv++) {fictiousVtx.push_back(fv); TS.grains[fv].id = fv; TS.grains[fv].isSphere = false;}}
+ if (fictiousVtx.size()==6){
CGT::Point& Pmin = TS.box.base; CGT::Point& Pmax = TS.box.sommet;
- Real FAR = 1e4;
TS.grains[fictiousVtx[0]].sphere = CGT::Sphere(CGT::Point(0.5*(Pmin.x()+Pmax.x()),Pmin.y()-FAR*(Pmax.x()-Pmin.x()),0.5*(Pmax.z()+Pmin.z())),FAR*(Pmax.x()-Pmin.x()));
TS.grains[fictiousVtx[1]].sphere = CGT::Sphere(CGT::Point(0.5*(Pmin.x()+Pmax.x()),Pmax.y()+FAR*(Pmax.x()-Pmin.x()),0.5*(Pmax.z()+Pmin.z())),FAR*(Pmax.x()-Pmin.x()));
TS.grains[fictiousVtx[2]].sphere =
@@ -154,7 +159,8 @@
CGT::Sphere(CGT::Point(0.5*(Pmin.x()+Pmax.x()),0.5*(Pmax.y()+Pmin.y()),Pmin.z()-FAR*(Pmax.y()-Pmin.y())),FAR*(Pmax.y()-Pmin.y()));
TS.grains[fictiousVtx[5]].sphere =
CGT::Sphere(CGT::Point(0.5*(Pmin.x()+Pmax.x()),0.5*(Pmax.y()+Pmin.y()),Pmax.z()+FAR*(Pmax.y()-Pmin.y())),FAR*(Pmax.y()-Pmin.y()));
- }
+ } else LOG_INFO(" the number of fictious vertices should be 0 or 6 usually");
+
InteractionContainer::iterator ii = scene->interactions->begin();
InteractionContainer::iterator iiEnd = scene->interactions->end();
for (; ii!=iiEnd ; ++ii) {
@@ -198,7 +204,7 @@
LOG_DEBUG("stress controller engine found");
triaxialCompressionEngine = YADE_PTR_CAST<TriaxialCompressionEngine> (*itFirst);}
}
- if (!triaxialCompressionEngine) LOG_ERROR("stress controller engine not found");}
+ if (!triaxialCompressionEngine) LOG_INFO("stress controller engine not found");}
if (triaxialCompressionEngine) {
TS.wszzh = triaxialCompressionEngine->stress[triaxialCompressionEngine->wall_top][1];
=== modified file 'pkg/dem/MicroMacroAnalyser.hpp'
--- pkg/dem/MicroMacroAnalyser.hpp 2010-11-07 11:46:20 +0000
+++ pkg/dem/MicroMacroAnalyser.hpp 2012-01-13 11:24:44 +0000
@@ -53,6 +53,7 @@
((int,interval,100,,"Number of timesteps between analyzed states."))
((bool,compDeformation,false,,"Is the engine just saving states or also computing and outputing deformations for each increment?"))
((bool,compIncrt,false,,"Should increments of force and displacements be defined on [n,n+1]? If not, states will be saved with only positions and forces (no displacements)."))
+ ((bool,nonSphereAsFictious,true,,"bodies that are not spheres will be used to defines bounds (else just skipped)."))
,/*init*/
,/*ctor*/
analyser = shared_ptr<CGT::KinematicLocalisationAnalyser> (new CGT::KinematicLocalisationAnalyser);
=== modified file 'pkg/dem/TesselationWrapper.cpp'
--- pkg/dem/TesselationWrapper.cpp 2012-01-11 21:30:18 +0000
+++ pkg/dem/TesselationWrapper.cpp 2012-01-13 11:24:44 +0000
@@ -117,12 +117,6 @@
//cerr << " loaded : " << Ng<<", triangulated : "<<TW.n_spheres<<", mean radius = " << TW.mean_radius<<endl;
}
-/*double pminx=0;
-double pminy=0;
-double pminz=0;
-double pmaxx=0;
-double pmaxy=0;
-double pmaxz=0;*/
double thickness = 0;
TesselationWrapper::~TesselationWrapper() { if (Tes) delete Tes;}
@@ -167,12 +161,7 @@
bool TesselationWrapper::insert(double x, double y, double z, double rad, unsigned int id)
{
using namespace std;
- Pmin = CGT::Point(min(Pmin.x(), x-rad),
- min(Pmin.y(), y-rad),
- min(Pmin.z(), z-rad));
- Pmax = CGT::Point(max(Pmax.x(), x+rad),
- max(Pmax.y(), y+rad),
- max(Pmax.z(), z+rad));
+ checkMinMax(x,y,z,rad);
mean_radius += rad;
++n_spheres;
return (Tes->insert(x,y,z,rad,id)!=NULL);
@@ -181,29 +170,15 @@
void TesselationWrapper::checkMinMax(double x, double y, double z, double rad)
{
using namespace std;
- Pmin = CGT::Point(min(Pmin.x(), x-rad),
- min(Pmin.y(), y-rad),
- min(Pmin.z(), z-rad));
- Pmax = CGT::Point(max(Pmax.x(), x+rad),
- max(Pmax.y(), y+rad),
- max(Pmax.z(), z+rad));
- mean_radius += rad;
- ++n_spheres;
+ Pmin = CGT::Point(min(Pmin.x(), x-rad), min(Pmin.y(), y-rad), min(Pmin.z(), z-rad));
+ Pmax = CGT::Point(max(Pmax.x(), x+rad), max(Pmax.y(), y+rad), max(Pmax.z(), z+rad));
}
bool TesselationWrapper::move(double x, double y, double z, double rad, unsigned int id)
{
using namespace std;
-
- Pmin = CGT::Point(min(Pmin.x(), x-rad),
- min(Pmin.y(), y-rad),
- min(Pmin.z(), z-rad));
- Pmax = CGT::Point(max(Pmax.x(), x+rad),
- max(Pmax.y(), y+rad),
- max(Pmax.z(), z+rad));
- mean_radius += rad;
-
+ checkMinMax(x,y,z,rad);
if (Tes->move(x,y,z,rad,id)!=NULL)
return true;
else {
=== modified file 'pkg/dem/TesselationWrapper.hpp'
--- pkg/dem/TesselationWrapper.hpp 2012-01-11 21:30:18 +0000
+++ pkg/dem/TesselationWrapper.hpp 2012-01-13 11:24:44 +0000
@@ -108,6 +108,7 @@
facet_end = Tes->Triangulation().finite_edges_end();
facet_it = Tes->Triangulation().finite_edges_begin();
inf=1e10;
+ mma.analyser->SetConsecutive(false);
,/*py*/
.def("triangulate",&TesselationWrapper::insertSceneSpheres,(python::arg("reset")=true),"triangulate spheres of the packing")
.def("setState",&TesselationWrapper::setState,(python::arg("state")=0),"Make the current state the initial (0) or final (1) configuration for the definition of displacement increments, use only state=0 if you just want to get only volmumes and porosity.")
_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to : yade-dev@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-dev
More help : https://help.launchpad.net/ListHelp
