yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #11551
[Branch ~yade-pkg/yade/git-trunk] Rev 3491: Added volumePower attribute to Law2_PolyhedraGeom_PolyhedraPhys_Volumetric, Polyhedra code small ...
------------------------------------------------------------
revno: 3491
committer: Jan Stransky <jan.stransky@xxxxxxxxxxx>
timestamp: Mon 2014-10-20 10:28:15 +0200
message:
Added volumePower attribute to Law2_PolyhedraGeom_PolyhedraPhys_Volumetric, Polyhedra code small cleanup
modified:
pkg/dem/Polyhedra.cpp
pkg/dem/Polyhedra.hpp
pkg/dem/Polyhedra_Ig2.cpp
pkg/dem/Polyhedra_support.cpp
--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== modified file 'pkg/dem/Polyhedra.cpp'
--- pkg/dem/Polyhedra.cpp 2014-10-17 15:57:30 +0000
+++ pkg/dem/Polyhedra.cpp 2014-10-20 08:28:15 +0000
@@ -97,7 +97,7 @@
std::transform( P.points_begin(), P.points_end(), P.points_begin(), t_trans);
//compute inertia
- double vtet;
+ Real vtet;
Vector3r ctet;
Matrix3r Itet1, Itet2;
Matrix3r inertia_tensor(Matrix3r::Zero());
@@ -126,8 +126,8 @@
// I_new = eigenvalues on diagonal
// set positove direction of vectors - otherwise reloading does not work
Matrix3r sign(Matrix3r::Zero());
- double max_v_signed = I_rot(0,0);
- double max_v = std::abs(I_rot(0,0));
+ Real max_v_signed = I_rot(0,0);
+ Real max_v = std::abs(I_rot(0,0));
if (max_v < std::abs(I_rot(1,0))) {max_v_signed = I_rot(1,0); max_v = std::abs(I_rot(1,0));}
if (max_v < std::abs(I_rot(2,0))) {max_v_signed = I_rot(2,0); max_v = std::abs(I_rot(2,0));}
sign(0,0) = max_v_signed/max_v;
@@ -175,11 +175,11 @@
int iter = 0;
bool isOK;
//fill box 5x5x5 with randomly located nuclei with restricted minimal mutual distance 0.75 which gives approximate mean mutual distance 1;
- double dist_min2 = 0.75*0.75;
+ Real dist_min2 = 0.75*0.75;
while(iter<500){
isOK = true;
iter++;
- trial = CGALpoint(double(rand())/RAND_MAX*5.+2.5,double(rand())/RAND_MAX*5.+2.5,double(rand())/RAND_MAX*5.+2.5);
+ trial = CGALpoint(Real(rand())/RAND_MAX*5.+2.5,Real(rand())/RAND_MAX*5.+2.5,Real(rand())/RAND_MAX*5.+2.5);
for(int i=0;i< (int) nuclei.size();i++) {
isOK = pow(nuclei[i].x()-trial.x(),2)+pow(nuclei[i].y()-trial.y(),2)+pow(nuclei[i].z()-trial.z(),2) > dist_min2;
if (!isOK) break;
@@ -205,7 +205,7 @@
//resize and rotate the voronoi cell
- Quaternionr Rot(double(rand())/RAND_MAX,double(rand())/RAND_MAX,double(rand())/RAND_MAX,double(rand())/RAND_MAX);
+ Quaternionr Rot(Real(rand())/RAND_MAX,Real(rand())/RAND_MAX,Real(rand())/RAND_MAX,Real(rand())/RAND_MAX);
Rot.normalize();
for(int i=0; i< (int) v.size();i++) {
v[i] = Rot*(Vector3r(v[i][0]*size[0],v[i][1]*size[1],v[i][2]*size[2]));
@@ -446,7 +446,8 @@
//zero penetration depth means no interaction force
if(!(contactGeom->equivalentPenetrationDepth > 1E-18) || !(contactGeom->penetrationVolume > 0)) {phys->normalForce = Vector3r(0.,0.,0.); phys->shearForce = Vector3r(0.,0.,0.); return true;}
- Vector3r normalForce=contactGeom->normal*contactGeom->penetrationVolume*phys->kn;
+ Real prop = std::pow(contactGeom->penetrationVolume,volumePower);
+ Vector3r normalForce=contactGeom->normal*prop*phys->kn;
//shear force: in case the polyhdras are separated and come to contact again, one should not use the previous shear force
if (contactGeom->isShearNew)
=== modified file 'pkg/dem/Polyhedra.hpp'
--- pkg/dem/Polyhedra.hpp 2014-10-17 15:57:30 +0000
+++ pkg/dem/Polyhedra.hpp 2014-10-20 08:28:15 +0000
@@ -58,7 +58,7 @@
void Initialize();
bool IsInitialized(){return init;}
std::vector<Vector3r> GetOriginalVertices();
- double GetVolume(){Initialize(); return volume;}
+ Real GetVolume(){Initialize(); return volume;}
Quaternionr GetOri(){Initialize(); return orientation;}
Polyhedron GetPolyhedron(){return P;};
void Clear(){v.clear(); P.clear(); init = 0; size = Vector3r(1.,1.,1.); faceTri.clear();};
@@ -73,7 +73,7 @@
//sign of performed initialization
bool init;
//centroid Volume
- double volume;
+ Real volume;
//centroid inerta - diagonal of the tensor
Vector3r inertia;
//orientation, that provides diagonal inertia tensor
@@ -149,15 +149,15 @@
/*! Elastic material */
class PolyhedraMat: public Material{
public:
- PolyhedraMat(double N, double S, double F){Kn=N; Ks=S; frictionAngle=F;};
- double GetStrength(){return strength;};
+ PolyhedraMat(Real N, Real S, Real F){Kn=N; Ks=S; frictionAngle=F;};
+ Real GetStrength(){return strength;};
virtual ~PolyhedraMat(){};
YADE_CLASS_BASE_DOC_ATTRS_CTOR(PolyhedraMat,Material,"Elastic material with Coulomb friction.",
- ((Real,Kn,1e8,,"Normal volumetric 'stiffness' (N/m3)."))
+ ((Real,Kn,1e8,,"Normal 'stiffness' (N/m3 for Law_.._Volumetric, N/m for Law2_.._Simple)."))
((Real,Ks,1e5,,"Shear stiffness (N/m)."))
((Real,frictionAngle,.5,,"Contact friction angle (in radians)."))
((bool,IsSplitable,0,,"To be splitted ... or not"))
- ((double,strength,100,,"Stress at whis polyhedra of volume 4/3*pi [mm] breaks.")),
+ ((Real,strength,100,,"Stress at whis polyhedra of volume 4/3*pi [mm] breaks.")),
/*ctor*/ createIndex();
);
REGISTER_CLASS_INDEX(PolyhedraMat,Material);
@@ -240,7 +240,7 @@
REGISTER_SERIALIZABLE(Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys);
//***************************************************************************
-/*! Calculate physical response based on penetration configuration given by TTetraGeom. */
+/*! Calculate physical response based on penetration configuration given by PolyhedraGeom. */
class Law2_PolyhedraGeom_PolyhedraPhys_Volumetric: public LawFunctor{
OpenMPAccumulator<Real> plasticDissipation;
@@ -248,14 +248,15 @@
Real elasticEnergy ();
Real getPlasticDissipation();
void initPlasticDissipation(Real initVal=0);
- YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_PolyhedraGeom_PolyhedraPhys_Volumetric,LawFunctor,"Calculate physical response of 2 :yref:`vector<Polyhedra>` in interaction, based on penetration configuration given by :yref:`PolyhedraGeom`.",
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_PolyhedraGeom_PolyhedraPhys_Volumetric,LawFunctor,"Calculate physical response of 2 :yref:`vector<Polyhedra>` in interaction, based on penetration configuration given by :yref:`PolyhedraGeom`. Normal force is proportional to the volume of intersection",
+ ((Real,volumePower,1.,,"Power of volume used in evaluation of normal force. Default is 1.0 - normal force is linearly proportional to volume. 1.0/3.0 would mean that normal force is proportional to the cube root of volume, approximation of penetration depth."))
((Vector3r,shearForce,Vector3r::Zero(),,"Shear force from last step"))
((bool,traceEnergy,false,,"Define the total energy dissipated in plastic slips at all contacts. This will trace only plastic energy in this law, see O.trackEnergy for a more complete energies tracing"))
((int,plastDissipIx,-1,(Attr::hidden|Attr::noSave),"Index for plastic dissipation (with O.trackEnergy)"))
((int,elastPotentialIx,-1,(Attr::hidden|Attr::noSave),"Index for elastic potential energy (with O.trackEnergy)"))
,,
.def("elasticEnergy",&Law2_PolyhedraGeom_PolyhedraPhys_Volumetric::elasticEnergy,"Compute and return the total elastic energy in all \"FrictPhys\" contacts")
- .def("plasticDissipation",&Law2_PolyhedraGeom_PolyhedraPhys_Volumetric::getPlasticDissipation,"Total energy dissipated in plastic slips at all FrictPhys contacts. Computed only if :yref:`Law2_ScGeom_FrictPhys_CundallStrack::traceEnergy` is true.")
+ .def("plasticDissipation",&Law2_PolyhedraGeom_PolyhedraPhys_Volumetric::getPlasticDissipation,"Total energy dissipated in plastic slips at all FrictPhys contacts. Computed only if :yref:`Law2_PolyhedraGeom_PolyhedraPhys_Volumetric::traceEnergy` is true.")
.def("initPlasticDissipation",&Law2_PolyhedraGeom_PolyhedraPhys_Volumetric::initPlasticDissipation,"Initialize cummulated plastic dissipation to a value (0 by default).")
);
FUNCTOR2D(PolyhedraGeom,PolyhedraPhys);
@@ -285,9 +286,9 @@
//Test if point is inside Polyhedron
bool Is_inside_Polyhedron(Polyhedron P, CGALpoint inside);
//return approximate intersection of sphere & polyhedron
-bool Sphere_Polyhedron_intersection(Polyhedron A, double r, CGALpoint C, CGALpoint centroid, double volume, CGALvector normal, double area);
+bool Sphere_Polyhedron_intersection(Polyhedron A, Real r, CGALpoint C, CGALpoint centroid, Real volume, CGALvector normal, Real area);
//return volume and centroid of polyhedra
-bool P_volume_centroid(Polyhedron P, double * volume, Vector3r * centroid);
+bool P_volume_centroid(Polyhedron P, Real * volume, Vector3r * centroid);
//CGAL - miniEigen communication
Vector3r FromCGALPoint(CGALpoint A);
Vector3r FromCGALVector(CGALvector A);
@@ -297,14 +298,14 @@
bool do_intersect(Polyhedron A, Polyhedron B);
bool do_intersect(Polyhedron A, Polyhedron B, std::vector<int> &sep_plane);
//connect triagular facets if possible
-Polyhedron Simplify(Polyhedron P, double lim);
+Polyhedron Simplify(Polyhedron P, Real lim);
//list of facets and edges
void PrintPolyhedron(Polyhedron P);
void PrintPolyhedron2File(Polyhedron P,FILE* X);
//normal by least square fitting of separating segments
Vector3r FindNormal(Polyhedron Int, Polyhedron PA, Polyhedron PB);
//calculate area of projection of polyhedron into the plane
-double CalculateProjectionArea(Polyhedron Int, CGALvector CGALnormal);
+Real CalculateProjectionArea(Polyhedron Int, CGALvector CGALnormal);
//split polyhedron
shared_ptr<Body> SplitPolyhedra(const shared_ptr<Body>& body, Vector3r direction, Vector3r point);
//new polyhedra
=== modified file 'pkg/dem/Polyhedra_Ig2.cpp'
--- pkg/dem/Polyhedra_Ig2.cpp 2014-09-22 12:28:50 +0000
+++ pkg/dem/Polyhedra_Ig2.cpp 2014-10-20 08:28:15 +0000
@@ -60,7 +60,7 @@
Int = Polyhedron_Polyhedron_intersection(PA,PB,ToCGALPoint(bang->contactPoint),ToCGALPoint(se31.position),ToCGALPoint(se32.position+shift2), bang->sep_plane);
//volume and centroid of intersection
- double volume;
+ Real volume;
Vector3r centroid;
P_volume_centroid(Int, &volume, ¢roid);
if(isnan(volume) || volume<=1E-25 || volume > min(A->GetVolume(),B->GetVolume())) {
@@ -76,10 +76,10 @@
if((se32.position+shift2-centroid).dot(normal)<0) normal*=-1;
//calculate area of projection of Intersection into the normal plane
- //double area = CalculateProjectionArea(Int, ToCGALVector(normal));
+ //Real area = CalculateProjectionArea(Int, ToCGALVector(normal));
//if(isnan(area) || area<=1E-20) {bang->equivalentPenetrationDepth=0; return true;}
- double area = volume/1E-8;
-
+ //Real area = volume/1E-8;
+ Real area = std::pow(volume,2./3.);
// store calculated stuff in bang; some is redundant
bang->equivalentCrossSection=area;
bang->contactPoint=centroid;
@@ -145,15 +145,15 @@
Int = Polyhedron_Plane_intersection(PB,A,ToCGALPoint(se32.position),ToCGALPoint(bang->contactPoint));
//volume and centroid of intersection
- double volume;
+ Real volume;
Vector3r centroid;
P_volume_centroid(Int, &volume, ¢roid);
if(isnan(volume) || volume<=1E-25 || volume > B->GetVolume()) {bang->equivalentPenetrationDepth=0; return true;}
if (!Is_inside_Polyhedron(PB, ToCGALPoint(centroid))) {bang->equivalentPenetrationDepth=0; return true;}
//calculate area of projection of Intersection into the normal plane
- double area = volume/1E-8;
- //double area = CalculateProjectionArea(Int, CGALnormal);
+ Real area = volume/1E-8;
+ //Real area = CalculateProjectionArea(Int, CGALnormal);
//if(isnan(area) || area<=1E-20) {bang->equivalentPenetrationDepth=0; return true;}
// store calculated stuff in bang; some is redundant
@@ -206,7 +206,7 @@
v[1]=help;
}
- double f_area = sqrt(f_normal.squared_length());
+ Real f_area = sqrt(f_normal.squared_length());
for (int i=3; i<6; i++) v[i] = v[i-3]-f_normal/f_area*0.05*sqrt(f_area); // vertices in global coordinates
Polyhedron PA;
@@ -240,7 +240,7 @@
Int = Polyhedron_Polyhedron_intersection(PA,PB,ToCGALPoint(bang->contactPoint),ToCGALPoint(se31.position),ToCGALPoint(se32.position), bang->sep_plane);
//volume and centroid of intersection
- double volume;
+ Real volume;
Vector3r centroid;
P_volume_centroid(Int, &volume, ¢roid);
if(isnan(volume) || volume<=1E-25 || volume > B->GetVolume()) {bang->equivalentPenetrationDepth=0; return true;}
@@ -251,8 +251,8 @@
if((se32.position-centroid).dot(normal)<0) normal*=-1;
//calculate area of projection of Intersection into the normal plane
- double area = volume/1E-8;
- //double area = CalculateProjectionArea(Int, ToCGALVector(normal));
+ Real area = volume/1E-8;
+ //Real area = CalculateProjectionArea(Int, ToCGALVector(normal));
//if(isnan(area) || area<=1E-20) {bang->equivalentPenetrationDepth=0; return true;}
// store calculated stuff in bang; some is redundant
@@ -301,7 +301,7 @@
}
//volume and centroid of intersection
- double volume, area;
+ Real volume, area;
CGALvector normalCGAL;
CGALpoint centroidCGAL=ToCGALPoint(se32.position);
Sphere_Polyhedron_intersection(PB, r, ToCGALPoint(se31.position), centroidCGAL, volume, normalCGAL, area);
=== modified file 'pkg/dem/Polyhedra_support.cpp'
--- pkg/dem/Polyhedra_support.cpp 2014-07-03 07:16:58 +0000
+++ pkg/dem/Polyhedra_support.cpp 2014-10-20 08:28:15 +0000
@@ -24,13 +24,13 @@
//**********************************************************************************
//return volume and centroid of polyhedron
-bool P_volume_centroid(Polyhedron P, double * volume, Vector3r * centroid){
+bool P_volume_centroid(Polyhedron P, Real * volume, Vector3r * centroid){
Vector3r basepoint = FromCGALPoint(P.vertices_begin()->point());
Vector3r A,B,C,D;
(*volume) = 0;
- double vtet;
+ Real vtet;
(*centroid) = Vector3r(0.,0.,0.);
//compute centroid and volume
@@ -77,26 +77,26 @@
#define z4 dv[2]
// Jacobian of transformation to the reference 4hedron
- double detJ=(x2-x1)*(y3-y1)*(z4-z1)+(x3-x1)*(y4-y1)*(z2-z1)+(x4-x1)*(y2-y1)*(z3-z1)
+ Real detJ=(x2-x1)*(y3-y1)*(z4-z1)+(x3-x1)*(y4-y1)*(z2-z1)+(x4-x1)*(y2-y1)*(z3-z1)
-(x2-x1)*(y4-y1)*(z3-z1)-(x3-x1)*(y2-y1)*(z4-z1)-(x4-x1)*(y3-y1)*(z2-z1);
detJ=std::abs(detJ);
- double a=detJ*(y1*y1+y1*y2+y2*y2+y1*y3+y2*y3+
+ Real a=detJ*(y1*y1+y1*y2+y2*y2+y1*y3+y2*y3+
y3*y3+y1*y4+y2*y4+y3*y4+y4*y4+z1*z1+z1*z2+
z2*z2+z1*z3+z2*z3+z3*z3+z1*z4+z2*z4+z3*z4+z4*z4)/60.;
- double b=detJ*(x1*x1+x1*x2+x2*x2+x1*x3+x2*x3+x3*x3+
+ Real b=detJ*(x1*x1+x1*x2+x2*x2+x1*x3+x2*x3+x3*x3+
x1*x4+x2*x4+x3*x4+x4*x4+z1*z1+z1*z2+z2*z2+z1*z3+
z2*z3+z3*z3+z1*z4+z2*z4+z3*z4+z4*z4)/60.;
- double c=detJ*(x1*x1+x1*x2+x2*x2+x1*x3+x2*x3+x3*x3+x1*x4+
+ Real c=detJ*(x1*x1+x1*x2+x2*x2+x1*x3+x2*x3+x3*x3+x1*x4+
x2*x4+x3*x4+x4*x4+y1*y1+y1*y2+y2*y2+y1*y3+
y2*y3+y3*y3+y1*y4+y2*y4+y3*y4+y4*y4)/60.;
// a' in the article etc.
- double a__=detJ*(2*y1*z1+y2*z1+y3*z1+y4*z1+y1*z2+
+ Real a__=detJ*(2*y1*z1+y2*z1+y3*z1+y4*z1+y1*z2+
2*y2*z2+y3*z2+y4*z2+y1*z3+y2*z3+2*y3*z3+
y4*z3+y1*z4+y2*z4+y3*z4+2*y4*z4)/120.;
- double b__=detJ*(2*x1*z1+x2*z1+x3*z1+x4*z1+x1*z2+
+ Real b__=detJ*(2*x1*z1+x2*z1+x3*z1+x4*z1+x1*z2+
2*x2*z2+x3*z2+x4*z2+x1*z3+x2*z3+2*x3*z3+
x4*z3+x1*z4+x2*z4+x3*z4+2*x4*z4)/120.;
- double c__=detJ*(2*x1*y1+x2*y1+x3*y1+x4*y1+x1*y2+
+ Real c__=detJ*(2*x1*y1+x2*y1+x3*y1+x4*y1+x1*y2+
2*x2*y2+x3*y2+x4*y2+x1*y3+x2*y3+2*x3*y3+
x4*y3+x1*y4+x2*y4+x3*y4+2*x4*y4)/120.;
@@ -122,8 +122,8 @@
//**********************************************************************************
//distace of point from a plane (squared) with sign
-double Oriented_squared_distance(Plane P, CGALpoint x){
- double h = P.a()*x.x()+P.b()*x.y()+P.c()*x.z()+P.d();
+Real Oriented_squared_distance(Plane P, CGALpoint x){
+ Real h = P.a()*x.x()+P.b()*x.y()+P.c()*x.z()+P.d();
return ((h>0.)-(h<0.))*pow(h,2)/(CGALvector(P.a(),P.b(),P.c())).squared_length();
}
@@ -139,7 +139,7 @@
//**********************************************************************************
// test if point is inside polyhedra not closer than lim to its boundary
-bool Is_inside_Polyhedron(Polyhedron P, CGALpoint inside, double lim){
+bool Is_inside_Polyhedron(Polyhedron P, CGALpoint inside, Real lim){
Polyhedron::Plane_iterator pi;
lim = pow(lim,2);
for(pi=P.planes_begin(); pi!=P.planes_end(); ++pi){
@@ -200,7 +200,7 @@
Plane X(eIter1->vertex()->point(),CGAL::cross_product((eIter1->vertex()->point() - eIter1->opposite()->vertex()->point()),(eIter2->vertex()->point() - eIter2->opposite()->vertex()->point())));
if (!X.has_on_positive_side(B.vertices_begin()->point())) X = X.opposite();
- double lim = pow(DISTANCE_LIMIT,2);
+ Real lim = pow(DISTANCE_LIMIT,2);
for (Polyhedron::Vertex_iterator vIter = A.vertices_begin(); vIter != A.vertices_end(); ++vIter){
if (Oriented_squared_distance(X, vIter->point())>lim) {found = false; break;};
}
@@ -234,7 +234,7 @@
//test all pairs of edges from A & B
Plane X;
CGALvector vA;
- double lim = pow(DISTANCE_LIMIT,2);
+ Real lim = pow(DISTANCE_LIMIT,2);
i = 0;
for (Polyhedron::Edge_iterator eIter1 = A.edges_begin(); eIter1 != A.edges_end(); ++eIter1, i++){
vA = eIter1->vertex()->point() - eIter1->opposite()->vertex()->point();
@@ -259,9 +259,9 @@
//**********************************************************************************
//norm of difference between two planes
-double PlaneDifference(const Plane &a, const Plane &b){
- double la = sqrt(pow(a.a(),2) + pow(a.b(),2) + pow(a.c(),2) + pow(a.d(),2));
- double lb = sqrt(pow(b.a(),2) + pow(b.b(),2) + pow(b.c(),2) + pow(b.d(),2));
+Real PlaneDifference(const Plane &a, const Plane &b){
+ Real la = sqrt(pow(a.a(),2) + pow(a.b(),2) + pow(a.c(),2) + pow(a.d(),2));
+ Real lb = sqrt(pow(b.a(),2) + pow(b.b(),2) + pow(b.c(),2) + pow(b.d(),2));
return pow(a.a()/la-b.a()/lb,2) + pow(a.b()/la-b.b()/lb,2) + pow(a.c()/la-b.c()/lb,2) + pow(a.d()/la-b.d()/lb,2);
//in case we do not care of the orientation
@@ -270,7 +270,7 @@
//**********************************************************************************
//connect triagular facets if possible
-Polyhedron Simplify(Polyhedron P, double limit){
+Polyhedron Simplify(Polyhedron P, Real limit){
bool elimination = true;
while(elimination){
elimination = false;
@@ -367,7 +367,7 @@
//solve system of 3x3 by Cramers rule
Vector3r SolveLinSys3x3(Matrix3r A, Vector3r y){
//only system 3x3 by Cramers rule
- double det = A(0,0)*A(1,1)*A(2,2)+A(0,1)*A(1,2)*A(2,0)+A(0,2)*A(1,0)*A(2,1)-A(0,2)*A(1,1)*A(2,0)-A(0,1)*A(1,0)*A(2,2)-A(0,0)*A(1,2)*A(2,1);
+ Real det = A(0,0)*A(1,1)*A(2,2)+A(0,1)*A(1,2)*A(2,0)+A(0,2)*A(1,0)*A(2,1)-A(0,2)*A(1,1)*A(2,0)-A(0,1)*A(1,0)*A(2,2)-A(0,0)*A(1,2)*A(2,1);
if (det == 0) {cerr << "error in linear solver" << endl; return Vector3r(0,0,0);}
return Vector3r((y(0)*A(1,1)*A(2,2)+A(0,1)*A(1,2)*y(2)+A(0,2)*y(1)*A(2,1)-A(0,2)*A(1,1)*y(2)-A(0,1)*y(1)*A(2,2)-y(0)*A(1,2)*A(2,1))/det, (A(0,0)*y(1)*A(2,2)+y(0)*A(1,2)*A(2,0)+A(0,2)*A(1,0)*y(2)-A(0,2)*y(1)*A(2,0)-y(0)*A(1,0)*A(2,2)-A(0,0)*A(1,2)*y(2))/det, (A(0,0)*A(1,1)*y(2)+A(0,1)*y(1)*A(2,0)+y(0)*A(1,0)*A(2,1)-y(0)*A(1,1)*A(2,0)-A(0,1)*A(1,0)*y(2)-A(0,0)*y(1)*A(2,1))/det);
}
@@ -395,10 +395,10 @@
std::transform( PA.facets_begin(), PA.facets_end(), PA.planes_begin(),Plane_equation());
std::transform( PB.facets_begin(), PB.facets_end(), PB.planes_begin(),Plane_equation());
vector<bool> from_A(Int.size_of_facets());
- vector<double> minsA(Int.size_of_facets());
- vector<double> minsB(Int.size_of_facets());
+ vector<Real> minsA(Int.size_of_facets());
+ vector<Real> minsB(Int.size_of_facets());
int i=0;
- double minA, minB, k;
+ Real minA, minB, k;
for(pi = Int.planes_begin(); pi!=Int.planes_end(); ++pi,i++){
minA = 1.;
minB = 1.;
@@ -454,13 +454,13 @@
//**********************************************************************************
//not used
-double CalculateProjectionArea(Polyhedron Int, CGALvector CGALnormal){
+Real CalculateProjectionArea(Polyhedron Int, CGALvector CGALnormal){
//calculate area of projection of Intersection into the normal plane
- double area = 0.;
- double abs_size;
+ Real area = 0.;
+ Real abs_size;
Polyhedron::Halfedge_around_facet_circulator hfc0;
CGALvector normal2;
- double norm2;
+ Real norm2;
std::transform( Int.facets_begin(), Int.facets_end(), Int.planes_begin(),Plane_equation());
for(Polyhedron::Facet_iterator fi = Int.facets_begin(); fi!= Int.facets_end(); ++fi){
assert(fi->is_triangle());
@@ -477,7 +477,7 @@
//**********************************************************************************
//prepare data for CGAL convex hull
-vector<Plane> MergePlanes(vector<Plane> planes1, vector<Plane> planes2, double limit){
+vector<Plane> MergePlanes(vector<Plane> planes1, vector<Plane> planes2, Real limit){
vector<Plane> P = planes1;
bool add;
for(vector<Plane>::iterator i = planes2.begin(); i!=planes2.end(); ++i){
@@ -503,7 +503,7 @@
vector<CGALpoint> dual_planes;
// test if do intersect, find some intersecting point
bool intersection_found = false;
- double lim = pow(DISTANCE_LIMIT,2);
+ Real lim = pow(DISTANCE_LIMIT,2);
std::transform( A.facets_begin(), A.facets_end(), A.planes_begin(),Plane_equation());
// test centroid of previous intersection
if(Is_inside_Polyhedron(A,X,DISTANCE_LIMIT) && Oriented_squared_distance(B,X)<=-lim) {
@@ -603,10 +603,10 @@
} else {
if (!do_intersect(A, B, sep_plane)) return Intersection;
//some intersection point
- double dist_S, dist_T;
- double lim2 = pow(DISTANCE_LIMIT,2);
+ Real dist_S, dist_T;
+ Real lim2 = pow(DISTANCE_LIMIT,2);
CGALvector d1;
- double factor = sqrt(DISTANCE_LIMIT * 1.5);
+ Real factor = sqrt(DISTANCE_LIMIT * 1.5);
//test vertices A - not needed, edges are enough
for (Polyhedron::Vertex_iterator vIter = A.vertices_begin(); vIter != A.vertices_end() && !intersection_found; vIter++){
d1 = centroidA-vIter->point();
@@ -687,15 +687,15 @@
//**********************************************************************************
/*
//returnes intersection of Sphere and polyhedron (possibly empty) - not finished
-bool Sphere_Polyhedron_intersection(Polyhedron A, double r, CGALpoint C, CGALpoint centroid, double volume, CGALvector normal, double area){
+bool Sphere_Polyhedron_intersection(Polyhedron A, Real r, CGALpoint C, CGALpoint centroid, Real volume, CGALvector normal, Real area){
volume = 0;
area = 0;
- double r2 = pow(r,2);
+ Real r2 = pow(r,2);
//check points
Polyhedron::Vertex_iterator closest_vert;
- double dist2;
- double min_dist2 = (C-centroid).squared_length();
+ Real dist2;
+ Real min_dist2 = (C-centroid).squared_length();
//test vertices
for(Polyhedron::Vertex_iterator vi=A.vertices_begin(); vi!=A.vertices_end(); ++vi){
dist2 = (C-vi->point()).squared_length();
@@ -711,8 +711,8 @@
v1 = v1/sqrt(v1.squared_length());
Polyhedron::Halfedge_around_vertex_circulator closest_edge;
int i =0;
- double cos_angle;
- double max_cos_angle=-1;
+ Real cos_angle;
+ Real max_cos_angle=-1;
CGALvector v2;
for(Polyhedron::Halfedge_around_vertex_circulator hfi=closest_vert->vertex_begin(); i<(int)closest_vert->degree(); ++hfi, i++){
v2 = hfi->opposite()->vertex()->point()-closest_vert->point();
@@ -721,7 +721,7 @@
if(cos_angle > max_cos_angle) closest_edge=hfi; max_cos_angle=cos_angle; //find the edge with minimum angle
}
if(max_cos_angle<=0) return false; //all edges goes with angle above or equal 90
- double k = (C - closest_vert->point())*v2; //find closest point on the most descending edge (it must be on edge, otherwise opposite point would be already closer)
+ Real k = (C - closest_vert->point())*v2; //find closest point on the most descending edge (it must be on edge, otherwise opposite point would be already closer)
CGALpoint closest_point = closest_vert->point() + k*v2;
if((C-closest_point).squared_length()<r2){ //test if we are already inside the sphere
@@ -743,7 +743,7 @@
else closest_plane = closest_edge->opposite()->facet()->plane();
closest_point = closest_plane.projection(C);
if((C-closest_point).squared_length()<r2){ //test if we are already inside the sphere
- double h = sqrt((C-closest_point).squared_length());
+ Real h = sqrt((C-closest_point).squared_length());
volume = 3.1415*pow(h,2)/3*(3*r-h);
return true;
}
@@ -837,7 +837,7 @@
for(Polyhedron::Vertex_iterator vi = S2.vertices_begin(); vi != S2.vertices_end(); vi++){ v2.push_back(FromCGALPoint(vi->point()));};
shared_ptr<Body> BP = NewPolyhedra(v2, body->material);
//BP->shape->color = body->shape->color;
- BP->shape->color = Vector3r(double(rand())/RAND_MAX,double(rand())/RAND_MAX,double(rand())/RAND_MAX);
+ BP->shape->color = Vector3r(Real(rand())/RAND_MAX,Real(rand())/RAND_MAX,Real(rand())/RAND_MAX);
scene->bodies->insert(BP);
//set proper state variables
BP->state->vel = OrigVel + OrigAngVel.cross(BP->state->pos-OrigPos);
