Ig2_PP_PP_ScGeom.hpp 4.35 KB
Newer Older
1 2 3 4 5
/*CWBoon 2015 */
/* C.W. Boon, G.T. Houlsby, S. Utili (2013).  A new contact detection algorithm for three-dimensional non-spherical particles.  Powder Technology, 248, pp 94-102. */
/* code for calling MOSEK was for ver 6.  Please uncomment if you have the licence */

#pragma once
6 7 8 9 10 11 12
#ifdef YADE_POTENTIAL_PARTICLES
#include <lib/serialization/Serializable.hpp>
#include <pkg/dem/PotentialParticle.hpp>
#include <pkg/common/Dispatching.hpp>
#include <pkg/common/Sphere.hpp>
#include <Python.h>
#include <Eigen/Core>
13 14
#include <stdio.h>

15 16
# ifdef YADE_MOSEK
#include <mosek.h>
17 18 19 20 21 22
#include <math.h>
#include <stdlib.h>
#include <string.h>
#endif


23
class Ig2_PP_PP_ScGeom: public IGeomFunctor {
24

25 26
#ifdef YADE_MOSEK
	protected:
27 28 29 30 31 32 33 34 35
		std::string myfile;
		std::string Key;
		MSKrescodee mosekTaskEnv;
		MSKenv_t    mosekEnv;
#endif


	public :
		virtual bool go(const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& se32, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
36 37 38



39
		Real evaluatePP(const shared_ptr<Shape>& cm1, const State& state1, const Vector3r newTrial);
40
		void getPtOnParticle2(const shared_ptr<Shape>& cm1, const State& state1, Vector3r previousPt, Vector3r normal, Vector3r& newlocalPoint);
41

42
		bool contactPtMosekF2(const shared_ptr<Shape>& cm1, const State& state1, const shared_ptr<Shape>& cm2, const State& state2, Vector3r &contactPt);
43

44
		bool customSolve(const shared_ptr<Shape>& cm1, const State& state1, const shared_ptr<Shape>& cm2, const State& state2, Vector3r &contactPt, bool warmstart);
45 46


47
		Vector3r getNormal(const shared_ptr<Shape>& cm1, const State& state1, const Vector3r newTrial);
48

49 50 51 52 53 54
		void BrentZeroSurf(const shared_ptr<Shape>& cm1, const State& state1, const Vector3r bracketA, const Vector3r bracketB, Vector3r& zero);





55
		/////////////////////////////////////////////////////////////////////////////////////////////////////
56

57 58


59 60 61
		YADE_CLASS_BASE_DOC_ATTRS_CTOR(Ig2_PP_PP_ScGeom,IGeomFunctor,"EXPERIMENTAL. IGeom functor for PotentialParticle - PotentialParticle pair",
			((Real, accuracyTol, pow(10,-7),, "accuracy desired, tolerance criteria for SOCP"))
			((Real,interactionDetectionFactor,1.0,,"bool to avoid granular ratcheting")),
62 63 64 65
			//((std::string,myfile,"./PotentialParticles"+"","string")),
			//timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas);
			//mosekTaskEnv = MSK_makeenv(&mosekEnv,NULL,NULL,NULL,NULL);
			//mosekTaskEnv = MSK_initenv(mosekEnv);
66
		);
67 68 69



70 71 72 73 74

		FUNCTOR2D(PotentialParticle,PotentialParticle);
		// needed for the dispatcher, even if it is symmetric
		DEFINE_FUNCTOR_ORDER_2D(PotentialParticle,PotentialParticle);
		DECLARE_LOGGER;
75 76 77 78 79 80 81 82 83 84

};

REGISTER_SERIALIZABLE(Ig2_PP_PP_ScGeom);


#ifdef __cplusplus
extern "C" {
#endif

85
	/* LAPACK LU */
86 87
	//int dgesv(int varNo, int varNo2, Real *H, int varNo3, int *pivot, Real* g, int varNo4, int info){
	extern void dgesv_(const int *N, const int *nrhs, Real *Hessian, const int *lda, int *ipiv, Real *gradient, const int *ldb, int *info);
88 89 90
	// int ans;
	// dgesv_(&varNo, &varNo2, H, &varNo3, pivot,g, &varNo4, &ans);
	// return ans;
91
	//}
92

93
	/* LAPACK Cholesky */
94
	extern void dpbsv_(const char *uplo, const int *n, const int *kd, const int *nrhs, Real *AB, const int *ldab, Real *B, const int *ldb, int *info);
95

96
	/* LAPACK QR */
97
	extern void dgels_(const char *Trans, const int *m, const int *n, const int *nrhs, Real *A, const int *lda, Real *B, const int *ldb, const Real *work, const int *lwork, int *info);
98

99 100

	/*BLAS */
101
	extern void dgemm_(const char *transA, const char *transB, const int *m, const int *n, const int *k, const Real *alpha, Real *A, const int *lda, Real *B, const int *ldb, const Real *beta, Real *C, const int *ldc);
102

103
	extern void dgemv_(const char *trans, const int *m, const int *n, const Real *alpha, Real *A, const int *lda, Real *x, const int *incx, const Real *beta, Real *y, const int *incy);
104

105
	extern void dcopy_(const int *N, Real *x, const int *incx, Real *y, const int *incy);
106

107
	extern Real ddot_(const int *N, Real *x, const int *incx, Real *y, const int *incy);
108

109
	extern void daxpy_(const int *N, const Real *da, Real *dx, const int *incx, Real *dy, const int *incy);
110

111
	extern void dscal_(const int *N, const Real *alpha, Real *x, const int *incx);
112 113


114
	void dsyev_(const char *jobz, const char *uplo, const int *N, Real *A, const int *lda, Real *W, Real *work, int *lwork, int *info);
115 116 117 118 119 120


#ifdef __cplusplus
};
#endif

121
#endif // YADE_POTENTIAL_PARTICLES