yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #02129
[Branch ~yade-dev/yade/trunk] Rev 1769: 1. Hopefully fix vtk libs in SConstruct (?)
------------------------------------------------------------
revno: 1769
committer: Václav Šmilauer <vaclav@flux>
branch nick: trunk
timestamp: Sun 2009-10-04 16:37:41 +0200
message:
1. Hopefully fix vtk libs in SConstruct (?)
2. Add compress parameter to VTKRecorder (doesn't seem to make any difference in the output files, however...?)
3. Add CPM damage to VTKRecorder
4. Move spiralProject to Shop, to be usable from within c++
5. Add a 2 functions to compute 2d and 3d stresses to yade.eudoxos
modified:
SConstruct
lib/smoothing/WeightedAverage2d.hpp
pkg/common/DataClass/InteractingGeometry/InteractingFacet.hpp
pkg/dem/Engine/StandAloneEngine/VTKRecorder.cpp
pkg/dem/Engine/StandAloneEngine/VTKRecorder.hpp
pkg/dem/meta/Shop.cpp
pkg/dem/meta/Shop.hpp
py/SConscript
py/_eudoxos.cpp
py/_utils.cpp
--
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 'SConstruct'
--- SConstruct 2009-09-17 18:20:44 +0000
+++ SConstruct 2009-10-04 14:37:41 +0000
@@ -330,7 +330,8 @@
# check "optional" libs
if 'vtk' in env['features']:
conf.env.Append(CPPPATH=env['VTKINCDIR'])
- ok=conf.CheckLibWithHeader('vtkHybrid','vtkInstantiator.h','c++','vtkInstantiator::New();',autoadd=1)
+ ok=conf.CheckLibWithHeader(['vtkCommon'],'vtkInstantiator.h','c++','vtkInstantiator::New();',autoadd=1)
+ env.Append(LIBS='vtkHybrid')
if not ok: featureNotOK('vtk')
if 'opengl' in env['features']:
ok=conf.CheckLibWithHeader('glut','GL/glut.h','c++','glutGetModifiers();',autoadd=1)
=== modified file 'lib/smoothing/WeightedAverage2d.hpp'
--- lib/smoothing/WeightedAverage2d.hpp 2009-06-15 16:34:09 +0000
+++ lib/smoothing/WeightedAverage2d.hpp 2009-10-04 14:37:41 +0000
@@ -117,9 +117,9 @@
*
*/
template<typename T, typename Tvalue>
-struct SymmGaussianKernelAverage: public WeightedAverage<T,Tvalue> {
+struct SymmGaussianDistributionAverage: public WeightedAverage<T,Tvalue> {
Real stDev, relThreshold;
- SymmGaussianKernelAverage(const shared_ptr<GridContainer<T> >& _grid, Real _stDev, Real _relThreshold=3): WeightedAverage<T,Tvalue>(_grid), stDev(_stDev), relThreshold(_relThreshold){}
+ SymmGaussianDistributionAverage(const shared_ptr<GridContainer<T> >& _grid, Real _stDev, Real _relThreshold=3): WeightedAverage<T,Tvalue>(_grid), stDev(_stDev), relThreshold(_relThreshold){}
virtual Real getWeight(const Vector2r& meanPt, const T& e){
Vector2r pos=getPosition(e);
Real rSq=pow(meanPt[0]-pos[0],2)+pow(meanPt[1]-pos[1],2);
@@ -129,15 +129,15 @@
virtual vector<Vector2i> filterCells(const Vector2r& refPt){return WeightedAverage<T,Tvalue>::grid->circleFilter(refPt,stDev*relThreshold);}
};
-/* Class for doing template specialization of gaussian kernel average on SGKA_Scalar2d and for testing */
+/* Class for doing template specialization of gaussian kernel average on SGDA_Scalar2d and for testing */
struct Scalar2d{
Vector2r pos;
Real val;
};
/* Final specialization; it only needs to define getValue and getPosition -- these functions contain knowledge about the element class itself */
-struct SGKA_Scalar2d: public SymmGaussianKernelAverage<Scalar2d,Real>{
- SGKA_Scalar2d(const shared_ptr<GridContainer<Scalar2d> >& _grid, Real _stDev, Real _relThreshold=3): SymmGaussianKernelAverage<Scalar2d,Real>(_grid,_stDev,_relThreshold){}
+struct SGDA_Scalar2d: public SymmGaussianDistributionAverage<Scalar2d,Real>{
+ SGDA_Scalar2d(const shared_ptr<GridContainer<Scalar2d> >& _grid, Real _stDev, Real _relThreshold=3): SymmGaussianDistributionAverage<Scalar2d,Real>(_grid,_stDev,_relThreshold){}
virtual Real getValue(const Scalar2d& dp){return (Real)dp.val;}
virtual Vector2r getPosition(const Scalar2d& dp){return dp.pos;}
};
@@ -155,13 +155,13 @@
//struct Scalar2d{Vector2r pos; Real val;};
Vector2r tuple2vec2r(const python::tuple& t){return Vector2r(python::extract<Real>(t[0])(),python::extract<Real>(t[1])());}
Vector2i tuple2vec2i(const python::tuple& t){return Vector2i(python::extract<int>(t[0])(),python::extract<int>(t[1])());}
- shared_ptr<SGKA_Scalar2d> sgka;
+ shared_ptr<SGDA_Scalar2d> sgka;
struct Poly2d{vector<Vector2r> vertices; bool inclusive;};
vector<Poly2d> clips;
public:
pyGaussAverage(python::tuple lo, python::tuple hi, python::tuple nCells, Real stDev){
shared_ptr<GridContainer<Scalar2d> > g(new GridContainer<Scalar2d>(tuple2vec2r(lo),tuple2vec2r(hi),tuple2vec2i(nCells)));
- sgka=shared_ptr<SGKA_Scalar2d>(new SGKA_Scalar2d(g,stDev));
+ sgka=shared_ptr<SGDA_Scalar2d>(new SGDA_Scalar2d(g,stDev));
}
bool pointInsidePolygon(const Vector2r&,const vector<Vector2r>&);
bool ptIsClipped(const Vector2r& pt){
=== modified file 'pkg/common/DataClass/InteractingGeometry/InteractingFacet.hpp'
--- pkg/common/DataClass/InteractingGeometry/InteractingFacet.hpp 2009-07-17 20:50:55 +0000
+++ pkg/common/DataClass/InteractingGeometry/InteractingFacet.hpp 2009-10-04 14:37:41 +0000
@@ -8,7 +8,7 @@
#pragma once
-#include <yade/core/InteractingGeometry.hpp>
+#include<yade/core/InteractingGeometry.hpp>
#include<yade/core/Body.hpp>
// define this to have topology information about facets enabled;
=== modified file 'pkg/dem/Engine/StandAloneEngine/VTKRecorder.cpp'
--- pkg/dem/Engine/StandAloneEngine/VTKRecorder.cpp 2009-09-17 17:12:22 +0000
+++ pkg/dem/Engine/StandAloneEngine/VTKRecorder.cpp 2009-10-04 14:37:41 +0000
@@ -7,12 +7,15 @@
#include<vtkFloatArray.h>
#include<vtkUnstructuredGrid.h>
#include<vtkXMLUnstructuredGridWriter.h>
+#include<vtkZLibDataCompressor.h>
//#include<vtkXMLMultiBlockDataWriter.h>
//#include<vtkMultiBlockDataSet.h>
#include<vtkTriangle.h>
#include<yade/core/MetaBody.hpp>
#include<yade/pkg-common/Sphere.hpp>
#include<yade/pkg-common/Facet.hpp>
+#include<yade/pkg-dem/ConcretePM.hpp>
+
YADE_PLUGIN((VTKRecorder));
YADE_REQUIRE_FEATURE(VTK)
@@ -22,6 +25,7 @@
{
/* we always want to save the first state as well */
initRun=true;
+ compress=false;
}
VTKRecorder::~VTKRecorder()
@@ -40,7 +44,8 @@
if(rec=="spheres") recActive[REC_SPHERES]=true;
else if(rec=="facets") recActive[REC_FACETS]=true;
else if(rec=="colors") recActive[REC_COLORS]=true;
- else LOG_ERROR("Unknown recorder named `"<<rec<<"' (supported are: spheres, facets). Ignored.");
+ else if(rec=="cpmDamage") recActive[REC_CPM_DAMAGE]=true;
+ else LOG_ERROR("Unknown recorder named `"<<rec<<"' (supported are: spheres, facets, colors, cpmDamage). Ignored.");
}
vtkSmartPointer<vtkPoints> spheresPos = vtkSmartPointer<vtkPoints>::New();
@@ -51,6 +56,9 @@
vtkSmartPointer<vtkFloatArray> spheresColors = vtkSmartPointer<vtkFloatArray>::New();
spheresColors->SetNumberOfComponents(3);
spheresColors->SetName("Colors");
+ vtkSmartPointer<vtkFloatArray> damage = vtkSmartPointer<vtkFloatArray>::New();
+ damage->SetNumberOfComponents(1);
+ damage->SetName("cpmDamage");
vtkSmartPointer<vtkPoints> facetsPos = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkCellArray> facetsCells = vtkSmartPointer<vtkCellArray>::New();
@@ -58,6 +66,7 @@
facetsColors->SetNumberOfComponents(3);
facetsColors->SetName("Colors");
+
FOREACH(const shared_ptr<Body>& b, *rootBody->bodies){
if (recActive[REC_SPHERES])
{
@@ -75,6 +84,9 @@
float c[3] = {color[0],color[1],color[2]};
spheresColors->InsertNextTupleValue(c);
}
+ if (recActive[REC_CPM_DAMAGE]) {
+ damage->InsertNextValue(YADE_PTR_CAST<CpmMat>(b->physicalParameters)->normDmg);
+ }
continue;
}
}
@@ -105,6 +117,9 @@
}
}
+ vtkSmartPointer<vtkDataCompressor> compressor;
+ if(compress) compressor=vtkSmartPointer<vtkZLibDataCompressor>::New();
+
if (recActive[REC_SPHERES])
{
vtkSmartPointer<vtkUnstructuredGrid> spheresUg = vtkSmartPointer<vtkUnstructuredGrid>::New();
@@ -112,7 +127,9 @@
spheresUg->SetCells(VTK_VERTEX, spheresCells);
spheresUg->GetPointData()->AddArray(radii);
if (recActive[REC_COLORS]) spheresUg->GetPointData()->AddArray(spheresColors);
+ if (recActive[REC_CPM_DAMAGE]) spheresUg->GetPointData()->AddArray(damage);
vtkSmartPointer<vtkXMLUnstructuredGridWriter> writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
+ if(compress) writer->SetCompressor(compressor);
string fn=fileName+"spheres."+lexical_cast<string>(rootBody->currentIteration)+".vtu";
writer->SetFileName(fn.c_str());
writer->SetInput(spheresUg);
@@ -125,6 +142,7 @@
facetsUg->SetCells(VTK_TRIANGLE, facetsCells);
if (recActive[REC_COLORS]) facetsUg->GetCellData()->AddArray(facetsColors);
vtkSmartPointer<vtkXMLUnstructuredGridWriter> writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
+ if(compress) writer->SetCompressor(compressor);
string fn=fileName+"facets."+lexical_cast<string>(rootBody->currentIteration)+".vtu";
writer->SetFileName(fn.c_str());
writer->SetInput(facetsUg);
=== modified file 'pkg/dem/Engine/StandAloneEngine/VTKRecorder.hpp'
--- pkg/dem/Engine/StandAloneEngine/VTKRecorder.hpp 2009-09-16 04:58:04 +0000
+++ pkg/dem/Engine/StandAloneEngine/VTKRecorder.hpp 2009-10-04 14:37:41 +0000
@@ -4,17 +4,18 @@
class VTKRecorder: public PeriodicEngine {
public:
- enum {REC_SPHERES=0,REC_FACETS,REC_COLORS,REC_SENTINEL};
+ enum {REC_SPHERES=0,REC_FACETS,REC_COLORS,REC_CPM_DAMAGE,REC_SENTINEL};
//! A stuff to record: spheres,facets,colors
vector<string> recorders;
string fileName;
+ bool compress;
VTKRecorder();
~VTKRecorder();
void init(MetaBody*);
virtual void action(MetaBody*);
private:
- REGISTER_ATTRIBUTES(PeriodicEngine,(recorders)(fileName));
+ REGISTER_ATTRIBUTES(PeriodicEngine,(recorders)(fileName)(compress));
REGISTER_CLASS_AND_BASE(VTKRecorder,PeriodicEngine);
DECLARE_LOGGER;
};
=== modified file 'pkg/dem/meta/Shop.cpp'
--- pkg/dem/meta/Shop.cpp 2009-08-10 10:32:23 +0000
+++ pkg/dem/meta/Shop.cpp 2009-10-04 14:37:41 +0000
@@ -5,6 +5,7 @@
#include<boost/filesystem/convenience.hpp>
#include<boost/tokenizer.hpp>
+#include<boost/tuple/tuple.hpp>
#include<yade/core/MetaBody.hpp>
#include<yade/core/Body.hpp>
@@ -494,7 +495,40 @@
}
return dt;
}
-
+/* Project 3d point into 2d using spiral projection along given axis;
+ * the returned tuple is
+ *
+ * (height relative to the spiral, distance from axis, theta )
+ *
+ * dH_dTheta is the inclination of the spiral (height increase per radian),
+ * theta0 is the angle for zero height (by given axis).
+ */
+boost::tuple<Real,Real,Real> Shop::spiralProject(const Vector3r& pt, Real dH_dTheta, int axis, Real periodStart, Real theta0){
+ int ax1=(axis+1)%3,ax2=(axis+2)%3;
+ Real r=sqrt(pow(pt[ax1],2)+pow(pt[ax2],2));
+ Real theta;
+ if(r>Mathr::ZERO_TOLERANCE){
+ theta=acos(pt[ax1]/r);
+ if(pt[ax2]<0) theta=Mathr::TWO_PI-theta;
+ }
+ else theta=0;
+ Real hRef=dH_dTheta*(theta-theta0);
+ long period;
+ if(isnan(periodStart)){
+ Real h=Shop::periodicWrap(pt[axis]-hRef,hRef-Mathr::PI*dH_dTheta,hRef+Mathr::PI*dH_dTheta,&period);
+ return boost::make_tuple(r,h,theta);
+ }
+ else{
+ // Real hPeriodStart=(periodStart-theta0)*dH_dTheta;
+ //TRVAR4(hPeriodStart,periodStart,theta0,theta);
+ //Real h=Shop::periodicWrap(pt[axis]-hRef,hPeriodStart,hPeriodStart+2*Mathr::PI*dH_dTheta,&period);
+ theta=Shop::periodicWrap(theta,periodStart,periodStart+2*Mathr::PI,&period);
+ Real h=pt[axis]-hRef+period*2*Mathr::PI*dH_dTheta;
+ //TRVAR3(pt[axis],pt[axis]-hRef,period);
+ //TRVAR2(h,theta);
+ return boost::make_tuple(r,h,theta);
+ }
+}
shared_ptr<Interaction> Shop::createExplicitInteraction(body_id_t id1, body_id_t id2){
InteractionGeometryMetaEngine* geomMeta=NULL;
=== modified file 'pkg/dem/meta/Shop.hpp'
--- pkg/dem/meta/Shop.hpp 2009-08-10 10:32:23 +0000
+++ pkg/dem/meta/Shop.hpp 2009-10-04 14:37:41 +0000
@@ -83,6 +83,9 @@
//! Estimate timestep based on P-wave propagation speed
static Real PWaveTimeStep(shared_ptr<MetaBody> rb=shared_ptr<MetaBody>());
+ //! return 2d coordinates of a 3d point within plane defined by rotation axis and inclination of spiral, wrapped to the 0th period
+ static boost::tuple<Real, Real, Real> spiralProject(const Vector3r& pt, Real dH_dTheta, int axis=2, Real periodStart=std::numeric_limits<Real>::quiet_NaN(), Real theta0=0);
+
//! Calculate inscribed circle center of trianlge
static Vector3r inscribedCircleCenter(const Vector3r& v0, const Vector3r& v1, const Vector3r& v2);
=== modified file 'py/SConscript'
--- py/SConscript 2009-09-05 10:11:06 +0000
+++ py/SConscript 2009-10-04 14:37:41 +0000
@@ -11,7 +11,7 @@
]),
env.SharedLibrary('log',['log.cpp'],SHLIBPREFIX=''),
env.SharedLibrary('_utils',['_utils.cpp'],SHLIBPREFIX='',LIBS=env['LIBS']+[
- linkPlugins(['Shop','ConcretePM']),
+ linkPlugins(['Shop','ConcretePM','InteractingFacet']),
]),
env.SharedLibrary('_packPredicates',['_packPredicates.cpp'],SHLIBPREFIX='',
# link to the symlink to the python module (created in lib/SConstruct; see explanation there)
=== modified file 'py/_eudoxos.cpp'
--- py/_eudoxos.cpp 2009-08-20 09:38:22 +0000
+++ py/_eudoxos.cpp 2009-10-04 14:37:41 +0000
@@ -1,6 +1,7 @@
#include<yade/pkg-dem/ConcretePM.hpp>
#include<boost/python.hpp>
#include<yade/extra/boost_python_len.hpp>
+#include<yade/pkg-dem/Shop.hpp>
using namespace boost::python;
using namespace std;
#ifdef YADE_LOG4CXX
@@ -86,10 +87,66 @@
}
BOOST_PYTHON_FUNCTION_OVERLOADS(velocityTowardsAxis_overloads,velocityTowardsAxis,3,5);
+/* Compute σxx,σyy,σxy stresses over all spheres, in plane passing through given axis,
+ which will be coincident with the y axis in the 2d projection.
+ Not sure how much is this function useful... */
+std::vector<Vector3r> spiralSphereStresses2d(Real dH_dTheta,const int axis=2){
+ MetaBody* rb=Omega::instance().getRootBody().get();
+ vector<Vector3r> ret(rb->bodies->size(),Vector3r::ZERO);
+ int ax1=(axis+1)%3,ax2=(axis+2)%3;
+ FOREACH(const shared_ptr<Interaction>& I, *rb->interactions){
+ if(!I->isReal()) continue;
+ CpmPhys* phys=YADE_CAST<CpmPhys*>(I->interactionPhysics.get());
+ Dem3DofGeom* geom=YADE_CAST<Dem3DofGeom*>(I->interactionGeometry.get());
+ // get force in this interaction, project it to the plane we need
+ Vector3r force=phys->normalForce+phys->shearForce;
+ Vector3r y2dIn3d(Vector3r::ZERO); y2dIn3d[axis]=1.;
+ Vector3r x2dIn3d(geom->contactPoint); x2dIn3d[axis]=0.; x2dIn3d.Normalize();
+ //Vector3r planeNormal=(x2dIn3d).Cross(y2dIn3d); planeNormal.Normalize(); force-=force.Dot(planeNormal)*planeNormal;
+ // get contact point in 2d
+ Vector2r C; Real theta;
+ boost::tie(C[0],C[1],theta)=Shop::spiralProject(geom->contactPoint,dH_dTheta,axis);
+ // get force in 2d (ff is already projected to the plane)
+ Vector2r _ff(force.Dot(x2dIn3d),force.Dot(y2dIn3d));
+ // get positions in 2d (height relative to C, as particle's position could wrap to other spiral period than the contact point, which we don't want here.
+ const Vector3r& aa(geom->se31.position); const Vector3r& bb(geom->se32.position);
+ Vector2r pos[]={Vector2r(sqrt(pow(aa[ax1],2)+pow(aa[ax2],2)),C[1]+(aa[axis]-geom->contactPoint[axis])),Vector2r(sqrt(pow(bb[ax1],2)+pow(bb[ax2],2)),C[1]+(bb[axis]-geom->contactPoint[axis]))};
+ Vector2r ff[]={_ff,-_ff}; body_id_t ids[]={I->getId1(),I->getId2()};
+ for(int i=0; i<2; i++){
+ // signs of tension/compression along the respective axes
+ //Vector2r sgn(pos[i][0]<C[0]?1.:-1.,pos[i][1]<C[1]?1.:-1.);
+ int sgn=(C-pos[i]).Dot(ff[i])>0?1:-1; // force in the same direction as vector away from particle: tension (positive)
+ ret[ids[i]][0]+=.5*sgn*abs(ff[i][0])/phys->crossSection; ret[ids[i]][1]+=.5*sgn*abs(ff[i][1])/phys->crossSection;
+ // divide by the length
+ Real torque=(pos[i]-C)[0]*ff[i][1]-(pos[i]-C)[1]*ff[i][0];
+ ret[ids[i]][2]+=.5*torque/phys->crossSection;
+ }
+ }
+ return ret;
+}
+
+std::vector<Real> particleConfinement(){
+ MetaBody* rb=Omega::instance().getRootBody().get();
+ vector<Real> ret(rb->bodies->size(),0.);
+ FOREACH(const shared_ptr<Interaction>& I, *rb->interactions){
+ if(!I->isReal()) continue;
+ CpmPhys* phys=YADE_CAST<CpmPhys*>(I->interactionPhysics.get());
+ Dem3DofGeom* geom=YADE_CAST<Dem3DofGeom*>(I->interactionGeometry.get());
+ Vector3r f[]={phys->normalForce,-phys->normalForce}; Vector3r pos[]={geom->se31.position,geom->se32.position}; body_id_t ids[]={I->getId1(),I->getId2()};
+ Real stress=phys->normalForce.Length()/phys->crossSection;
+ for(int i=0; i<2; i++){
+ int sgn=(geom->contactPoint-pos[i]).Dot(f[i])>0?1:-1;
+ ret[ids[i]]+=sgn*stress;
+ }
+ }
+ return ret;
+}
BOOST_PYTHON_MODULE(_eudoxos){
def("velocityTowardsAxis",velocityTowardsAxis,velocityTowardsAxis_overloads(args("axisPoint","axisDirection","timeToAxis","subtractDist","perturbation")));
def("yieldSigmaTMagnitude",yieldSigmaTMagnitude);
+ def("spiralSphereStresses2d",spiralSphereStresses2d,(python::arg("dH_dTheta"),python::arg("axis")=2));
+ def("particleConfinement",particleConfinement);
}
=== modified file 'py/_utils.cpp'
--- py/_utils.cpp 2009-08-20 19:55:52 +0000
+++ py/_utils.cpp 2009-10-04 14:37:41 +0000
@@ -7,6 +7,7 @@
#include<yade/pkg-dem/SpheresContactGeometry.hpp>
#include<yade/pkg-dem/DemXDofGeom.hpp>
#include<yade/pkg-dem/SimpleViscoelasticBodyParameters.hpp>
+#include<yade/pkg-common/InteractingFacet.hpp>
#include<yade/pkg-common/NormalShearInteractions.hpp>
#include<yade/lib-computational-geometry/Hull2d.hpp>
#include<cmath>
@@ -214,9 +215,9 @@
}
return ret;
}
-/* Sum forces actiong on bodies within mask.
+/* Sum forces acting on bodies within mask.
*
- * @param mask groupMask of matching bodies
+ * @param ids list of ids
* @param direction direction in which forces are summed
*
*/
@@ -233,6 +234,17 @@
return ret;
}
+/* Sum force acting on facets given by their ids in the sense of their respective normals. */
+Real sumFacetNormalForces(vector<body_id_t> ids){
+ shared_ptr<MetaBody> rb=Omega::instance().getRootBody(); rb->bex.sync();
+ Real ret=0;
+ FOREACH(const body_id_t id, ids){
+ InteractingFacet* f=YADE_CAST<InteractingFacet*>(Body::byId(id,rb)->interactingGeometry.get());
+ ret+=rb->bex.getForce(id).Dot(f->nf);
+ }
+ return ret;
+}
+
/* Set wire display of all/some/none bodies depending on the filter. */
void wireSome(string filter){
enum{none,all,noSpheres,unknown};
@@ -354,41 +366,12 @@
}
-/* Project 3d point into 2d using spiral projection along given axis;
- * the returned tuple is
- *
- * ( (height relative to the spiral, distance from axis), theta )
- *
- * dH_dTheta is the inclination of the spiral (height increase per radian),
- * theta0 is the angle for zero height (by given axis).
- */
python::tuple spiralProject(const Vector3r& pt, Real dH_dTheta, int axis=2, Real periodStart=std::numeric_limits<Real>::quiet_NaN(), Real theta0=0){
- int ax1=(axis+1)%3,ax2=(axis+2)%3;
- Real r=sqrt(pow(pt[ax1],2)+pow(pt[ax2],2));
- Real theta;
- if(r>Mathr::ZERO_TOLERANCE){
- theta=acos(pt[ax1]/r);
- if(pt[ax2]<0) theta=Mathr::TWO_PI-theta;
- }
- else theta=0;
- Real hRef=dH_dTheta*(theta-theta0);
- long period;
- if(isnan(periodStart)){
- Real h=Shop::periodicWrap(pt[axis]-hRef,hRef-Mathr::PI*dH_dTheta,hRef+Mathr::PI*dH_dTheta,&period);
- return python::make_tuple(python::make_tuple(r,h),theta);
- }
- else{
- // Real hPeriodStart=(periodStart-theta0)*dH_dTheta;
- //TRVAR4(hPeriodStart,periodStart,theta0,theta);
- //Real h=Shop::periodicWrap(pt[axis]-hRef,hPeriodStart,hPeriodStart+2*Mathr::PI*dH_dTheta,&period);
- theta=Shop::periodicWrap(theta,periodStart,periodStart+2*Mathr::PI,&period);
- Real h=pt[axis]-hRef+period*2*Mathr::PI*dH_dTheta;
- //TRVAR3(pt[axis],pt[axis]-hRef,period);
- //TRVAR2(h,theta);
- return python::make_tuple(python::make_tuple(r,h),theta);
- }
+ Real r,h,theta;
+ boost::tie(r,h,theta)=Shop::spiralProject(pt,dH_dTheta,axis,periodStart,theta0);
+ return python::make_tuple(python::make_tuple(r,h),theta);
}
-BOOST_PYTHON_FUNCTION_OVERLOADS(spiralProject_overloads,spiralProject,2,5);
+//BOOST_PYTHON_FUNCTION_OVERLOADS(spiralProject_overloads,spiralProject,2,5);
// for now, don't return anything, since we would have to include the whole yadeControl.cpp because of pyInteraction
void Shop__createExplicitInteraction(body_id_t id1, body_id_t id2){ (void) Shop::createExplicitInteraction(id1,id2);}
@@ -419,11 +402,12 @@
def("kineticEnergy",Shop__kineticEnergy);
def("sumBexForces",sumBexForces);
def("sumBexTorques",sumBexTorques);
+ def("sumFacetNormalForces",sumFacetNormalForces);
def("forcesOnPlane",forcesOnPlane);
def("forcesOnCoordPlane",forcesOnCoordPlane);
def("totalForceInVolume",Shop__totalForceInVolume,"Return summed forces on all interactions and average isotropic stiffness, as tuple (Vector3,float)");
def("createInteraction",Shop__createExplicitInteraction);
- def("spiralProject",spiralProject,spiralProject_overloads(args("axis","periodStart","theta0")));
+ def("spiralProject",spiralProject,(python::arg("pt"),python::arg("dH_dTheta"),python::arg("axis")=2,python::arg("periodStart")=std::numeric_limits<Real>::quiet_NaN(),python::arg("theta0")=0));
def("pointInsidePolygon",pointInsidePolygon);
def("scalarOnColorScale",Shop::scalarOnColorScale);
def("highlightNone",highlightNone);
