yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #06524
[Branch ~yade-dev/yade/trunk] Rev 2609: 1. Add documentation for Ig2_Sphere_Sphere_L3Geom_Inc
------------------------------------------------------------
revno: 2609
committer: Václav Šmilauer <eu@xxxxxxxx>
branch nick: yade
timestamp: Sun 2010-12-12 00:01:30 +0100
message:
1. Add documentation for Ig2_Sphere_Sphere_L3Geom_Inc
2. Add regression tests for PBC & L3Geom
3. Fix syntax in other docs as necessary
modified:
pkg/dem/L3Geom.cpp
pkg/dem/L3Geom.hpp
pkg/dem/NozzleFactory.hpp
pkg/dem/WirePM.hpp
py/_extraDocs.py
py/tests/pbc.py
--
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 'pkg/dem/L3Geom.cpp'
--- pkg/dem/L3Geom.cpp 2010-12-07 16:18:48 +0000
+++ pkg/dem/L3Geom.cpp 2010-12-11 23:01:30 +0000
@@ -104,9 +104,11 @@
// noRatch: take radius or current distance as the branch vector; see discussion in ScGeom::precompute (avoidGranularRatcheting)
Vector3r c1x=(noRatch ? ( r1*normal).eval() : (contPt-state1.pos).eval());
Vector3r c2x=(noRatch ? (-r2*normal).eval() : (contPt-state2.pos+shift2).eval());
+ //Vector3r state2velCorrected=state2.vel+(scene->isPeriodic?scene->cell->intrShiftVel(I->cellDist):Vector3r::Zero()); // velocity of the second particle, corrected with meanfield velocity if necessary
+ //cerr<<"correction "<<(scene->isPeriodic?scene->cell->intrShiftVel(I->cellDist):Vector3r::Zero())<<endl;
Vector3r relShearVel=(state2.vel+state2.angVel.cross(c2x))-(state1.vel+state1.angVel.cross(c1x));
// account for relative velocity of particles in different cell periods
- if(scene->isPeriodic) relShearVel+=scene->cell->velGrad*scene->cell->Hsize*I->cellDist.cast<Real>();
+ if(scene->isPeriodic) relShearVel+=scene->cell->intrShiftVel(I->cellDist);
// separate the tangential component (this is likely to be not stricly necessary, since u[0] is set directly later)
relShearVel-=avgNormal.dot(relShearVel)*avgNormal;
Vector3r relShearDu=relShearVel*scene->dt;
@@ -202,7 +204,8 @@
Vector3r c1x=contPt-state1.pos; // is not aligned with normal; Wall could rotate in a plate-like manner
Vector3r c2x=(noRatch ? (-radius*normal).eval() : (contPt-state2.pos+shift2).eval());
- Vector3r relShearVel=(state2.vel+state2.angVel.cross(c2x))-(state1.vel+state1.angVel.cross(c1x));
+ Vector3r state2velCorrected=state2.vel+(scene->isPeriodic?scene->cell->intrShiftVel(I->cellDist):Vector3r::Zero()); // velocity of the second particle, corrected with meanfield velocity if necessary
+ Vector3r relShearVel=(state2velCorrected+state2.angVel.cross(c2x))-(state1.vel+state1.angVel.cross(c1x));
// perhaps not necessary, see in Ig2_Sphere_Sphere_L3Geom_Inc
relShearVel-=normal.dot(relShearVel)*normal;
Vector3r relShearDu=relShearVel*scene->dt;
=== modified file 'pkg/dem/L3Geom.hpp'
--- pkg/dem/L3Geom.hpp 2010-12-07 13:38:35 +0000
+++ pkg/dem/L3Geom.hpp 2010-12-11 23:01:30 +0000
@@ -90,7 +90,15 @@
YADE_CLASS_BASE_DOC_ATTRS(Ig2_Sphere_Sphere_L3Geom_Inc,IGeomFunctor,"Incrementally compute :yref:`L3Geom` for contact of 2 spheres.\n\n.. note:: The initial value of *u[0]* (normal displacement) might be non-zero, with or without *distFactor*, since it is given purely by sphere's geometry. If you want to set \"equilibrium distance\", do it in the contact law as explained in :yref:`L3Geom.u0`.",
((bool,noRatch,true,,"See :yref:`ScGeom.avoidGranularRatcheting`."))
((Real,distFactor,1,,"Create interaction if spheres are not futher than distFactor*(r1+r2)."))
- ((int,approxMask,0,,"Selectively enable geometrical approximations (bitmask); add the values for approximations to be enabled.\n\n1: do not renormalize transformation matrix at every step\n2: use previous transformation to transform velocities (which are known at mid-steps), instead of mid-step transformation computed as quaternion slerp at t=0.5.\n4: do not take average (mid-step) normal when computing relative shear displacement, use previous value instead\n8: do not re-normalize average (mid-step) normal, if used.…\nBy default, the mask is zero and neither of these approximations is used."))
+ ((int,approxMask,0,,"Selectively enable geometrical approximations (bitmask); add the values for approximations to be enabled.\n\n"
+ "== ===============================================================\n"
+ "1 do not renormalize transformation matrix at every step\n"
+ "2 use previous transformation to transform velocities (which are known at mid-steps), instead of mid-step transformation computed as quaternion slerp at t=0.5.\n"
+ "4 do not take average (mid-step) normal when computing relative shear displacement, use previous value instead\n"
+ "8 do not re-normalize average (mid-step) normal, if used.…\n"
+ "== ===============================================================\n\n"
+ "By default, the mask is zero, wherefore none of these approximations is used.\n"
+ ))
);
FUNCTOR2D(Sphere,Sphere);
DEFINE_FUNCTOR_ORDER_2D(Sphere,Sphere);
=== modified file 'pkg/dem/NozzleFactory.hpp'
--- pkg/dem/NozzleFactory.hpp 2010-11-07 11:46:20 +0000
+++ pkg/dem/NozzleFactory.hpp 2010-12-11 23:01:30 +0000
@@ -31,7 +31,7 @@
class DragForceApplier: public GlobalEngine{
public: virtual void action();
- YADE_CLASS_BASE_DOC_ATTRS(DragForceApplier,GlobalEngine,"Apply `drag force <http://en.wikipedia.org/wiki/Drag_equation>`__ on particles, decelerating them proportionally to their linear velocities. The applied force reads\n\n.. math:: F_{d}=-\\frac{\\vec{v}}{|\\vec{v}|}\\frac{1}{2}\\rho|\\vec{v}|^2 C_d A\n\nwhere $\\rho$ is the medium density (:yref:`density<DragForceApplier.density>`), $v$ is particle's velocity, $A$ is particle projected area (disc), $C_d$ is the drag coefficient (0.47 for :yref:`Sphere`), \n\n.. note:: Drag force is only applied to spherical particles.\n\n.. warn:: Not tested.",
+ YADE_CLASS_BASE_DOC_ATTRS(DragForceApplier,GlobalEngine,"Apply `drag force <http://en.wikipedia.org/wiki/Drag_equation>`__ on particles, decelerating them proportionally to their linear velocities. The applied force reads\n\n.. math:: F_{d}=-\\frac{\\vec{v}}{|\\vec{v}|}\\frac{1}{2}\\rho|\\vec{v}|^2 C_d A\n\nwhere $\\rho$ is the medium density (:yref:`density<DragForceApplier.density>`), $v$ is particle's velocity, $A$ is particle projected area (disc), $C_d$ is the drag coefficient (0.47 for :yref:`Sphere`), \n\n.. note:: Drag force is only applied to spherical particles.\n\n.. warning:: Not tested.",
((Real,density,0,,"Density of the medium."))
);
};
=== modified file 'pkg/dem/WirePM.hpp'
--- pkg/dem/WirePM.hpp 2010-12-10 02:42:10 +0000
+++ pkg/dem/WirePM.hpp 2010-12-11 23:01:30 +0000
@@ -47,7 +47,7 @@
((vector<Vector2r>,strainStressValues,,Attr::triggerPostLoad,"Piecewise linear definition of the stress-strain curve by set of points (strain[-]>0,stress[Pa]>0) for one single wire. Tension only is considered and the point (0,0) is not needed!"))
((bool,isDoubleTwist,false,,"Type of the mesh. If true two particles of the same material which body ids differ by one will be considered as double-twisted interaction."))
((Real,lambdaEps,0.4,,"Parameter between 0 and 1 to reduce the failure strain of the double-twisted wire (as used by [Bertrand2008]_). [-]"))
- ((Real,lambdak,0.21,,"Parameter between 0 and 1 to compute the elastic stiffness of the double-twisted wire (as used by [Bertrand2008]_): $\\k^D=2(\\lambda_k k_h + (1-\\lambda_k)k^S)$. [-]"))
+ ((Real,lambdak,0.21,,"Parameter between 0 and 1 to compute the elastic stiffness of the double-twisted wire (as used by [Bertrand2008]_): $k^D=2(\\lambda_k k_h + (1-\\lambda_k)k^S)$. [-]"))
((Real,as,0,Attr::readonly,"Cross-section area of a single wire used for the computation of the limit normal contact forces. [m²]"))
,
createIndex();
=== modified file 'py/_extraDocs.py'
--- py/_extraDocs.py 2010-10-17 01:06:13 +0000
+++ py/_extraDocs.py 2010-12-11 23:01:30 +0000
@@ -108,80 +108,75 @@
Examples :ysrc:`scripts/test/peri3dController_example1` and :ysrc:`scripts/test/peri3dController_triaxialCompression` explain usage and inputs of Peri3dController, :ysrc:`scripts/test/peri3dController_shear` is an example of using shear components and also simulation on rotatd cell.
'''
-#wrapper.Peri3dController.__doc__=r'''
-#Class for controlling independently strain and stress tensor components of the periodic cell. :yref:`goal<Peri3dController.goal>` are the goal values, while :yref:`stressMask<Peri3dController.stressMask>` determines which components prescribe stress and which prescribe strain.
-#
-#Current stiffness and stress tensors are update at every step using formulas presented in [Kuhl2001]_:
-#
-#.. math::
- #:nowrap:
-#
- #\begin{align*}
- #\tens{\sig}&=\frac{1}{V}\sum_c l^c\left[\tens{N}^c F_N^c+{\tens{T}^c}^T \otimes F_T^c \right] & (34) \\
- #\tens{K}&=\frac{1}{V}\sum_c {l^c}^2\left[K_N^c \tens{N}^c\otimes\tens{N}^c+K_T^c{\tens{T}^c}^T\otimes\tens{T}^c\right] & (35)
- #\end{align*}
-#
-#where $c$ traverses all interactions within the :yref:`Cell` with volume $V$, $l$ the current distance between particle centers within the interaction, $F_N$ the normal force magnitude, $\vec{F_T}$ the shear force (in global coordinates), $K_N$ and $K_T$ are (current) normal and shear contact stiffnesses [N/m].
-#
-#The projection tensors $\tens{N}$ and $\tens{T}$ are ($n$ is the interaction unit normal, $\delta$ is the Kronecker delta)
-#
-#.. math::
- #:nowrap:
-#
- #\begin{align*}
- #\tens{N}&=\vec{n}\otimes\vec{n} \\
- #\tens{N}\otimes\tens{N}&=\vec{n}\otimes\vec{n}\otimes\vec{n}\otimes\vec{n}=n_i n_j n_k n_l \\
- #\tens{T}&=\vec{n}\otimes\tens{I}_{\rm sym}-\vec{n}\otimes\vec{n}\otimes\vec{n}\\
- #T_{ijk}&=n_l\frac{1}{2}(\delta_{lj}\delta_{ik}+\delta_{lk}\delta_{ij})-n_i n_j n_k=n_j\delta_{ik}/2+n_k\delta_{ij}/2-n_i n_j n_k \\
- #\tens{T}^T&=\tens{I}_{\rm sym}\otimes\vec{n}-\vec{n}\otimes\vec{n}\otimes\vec{n}\\
- #T_{ijk}^T&=\frac{1}{2}(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})n_l-n_i n_j n_k=\delta_{ik}n_j/2+n_i\delta_{jk}/2-n_i n_j n_k\\
- #\tens{T}^T\otimes\tens{T}=T^T_{ijm}T_{mkl}&=(\delta_{im}n_j/2+n_i\delta_{jm}/2-n_i n_j n_m)(n_k\delta_{ml}/2+n_l\delta_{mk}/2-n_m n_k n_l)= \\
- #&=n_jn_k\delta_{im}\delta_{ml}/4+n_jn_l\delta_{im}\delta_{mk}/4-n_jn_kn_ln_m\delta_{im}/2+ \\
- #&\quad+n_in_k\delta_{jm}\delta_{ml}/4+n_in_l\delta_{jm}\delta_{mk}/4-n_in_kn_ln_m\delta_{jm}/2-\\
- #&\quad-n_jn_kn_ln_m\delta_{im}/2-n_in_kn_ln_m\delta_{jm}/2+n_in_jn_kn_ln_mn_m= \\
- #&=n_jn_k\delta_{il}/4+n_jn_l\delta_{ik}/4+n_in_k\delta_{jl}/4+n_in_l\delta_{jk}/4- \\
- #&\quad-n_in_jn_kn_l/2-n_in_jn_kn_l/2-n_in_jn_kn_l/2-n_in_jn_kn_l/2+n_in_jn_kn_ln_mn_m= \\
- #&=n_jn_k\delta_{il}/4+n_jn_l\delta_{ik}/4+n_in_k\delta_{jl}/4+n_in_l\delta_{jk}/4-2n_in_jn_kn_l+n_in_jn_kn_ln_mn_m \\
- #&=n_jn_k\delta_{il}/4+n_jn_l\delta_{ik}/4+n_in_k\delta_{jl}/4+n_in_l\delta_{jk}/4-n_in_jn_kn_l \quad\hbox{(using \$n_m n_m=1\$)} %%% \$ must be escaped, since it would be dumbly replaced by :math:`...`
- #\end{align*}
-#
-#The global elastic stiffness tensor then reads (omitting the $c$ indices of per-contact symbols in the summation for readability)
-#
-#.. math::
- #
- #K_{ijkl}=\frac{1}{V}\sum_{c}{l}^2\left[K_N n_i n_j n_k n_l+K_T\left(\frac{n_j n_k \delta_{il}+n_j n_l\delta_{ik}+n_i n_k\delta_{jl}+n_i n_l\delta_{jk}}{4}-n_i n_j n_kn_l\right)\right].
-#
-#Global strain tensor $\tens{\eps}$ is symmetric and has 6 degrees of freedom, whereas the current transformation matrix $\mathcal{T}$ of the cell (:yref:`Cell.trsf`) has 9 independent components -- it has 3 additional for arbitrary rotation, which induces no strain. We perform `polar decomposition <http://en.wikipedia.org/wiki/Polar_decomposition#Matrix_polar_decomposition>`_ $\mathcal{T}=\mat{U}\mat{P}$, where $\mat{U}$ is unitary (representing rotation) and $\mat{P}$ is deformation in the proper sense; after converting $\mat{P}$ to strain, it is rotated back by $\mat{U}$ so that it is in global coordinates. To have `true strain <http://en.wikipedia.org/wiki/True_strain#True_strain>`_, the strain matrix is obtained via
-#
-#.. math:: \mat{\eps}=\mat{U}\log \mat{P}
-#
-#computed (see `wikipedia <http://en.wikipedia.org/wiki/Logarithm_of_a_matrix#Calculating_the_logarithm_of_a_diagonalizable_matrix>`_ and [Lu1998]_) by spectral decomposition
-#
-#.. math::
- #:nowrap:
- #
- #\begin{align*}
- #\mat{\eps}'&=\mat{V}^{-1}\mat{\eps}\mat{V} \\
- #\log(\mat{\eps})&=\mat{V}\log(\mat{\eps}')\mat{V}^{-1}
- #\end{align*}
-#
-#where $\mat{\eps}'$ is diagonal and $\log(\mat{\eps}')$ is a diagonal matrix with $\log(\mat{\eps}')_{ii}=\log(\mat{\eps}'_{ii})$. $\mat{V}$ is orthonormal, therefore $\mat{V}^{-1}=\mat{V}^T$.
-#
-#Once current values of $\tens{\sigma}$, $\tens{K}$ and $\tens{\eps}$ are known, the control algorithm comes. These tensors are converted to 6-vectors using the `Voigt notation <http://en.wikiversity.org/wiki/Introduction_to_Elasticity/Constitutive_relations#Voigt_notation>`_. Since for each component $i\in\{0\dots5\}$ either $\vec{\sig}_i$ of $\vec{\eps}_i$ is prescribed (determined by :yref:`stressMask<Peri3dController.stressMask>`), the equation $\vec{\sig}=\mat{K}\vec{\eps}$; assuming large strains, we take the rate form $\vec{\dot\sig}=\mat{K}\vec{\dot\eps}$ is decomposed in 2 parts:
-#
-#* $p$ denotes the part where strain component $\vec{\dot\eps}_i$ is *prescribed* (via $\frac{\vec{\eps}^g_i-\vec{\eps}_i}{\Dt}$, where $\vec{\eps}^g_i$ is the corresponding :yref:`goal<Peri3dController.goal>` strain component value);
-#* $u$ denotes the part where strain component $\vec{\dot\eps}_i$ is *unprescribed* (and $\vec{\sig}_i$ is prescribed, again via $\frac{\vec{\sig}^g_i-\vec{\sig}_i}{\Dt}$).
-#
-#The equation then is re-arranged as
-#
-#.. math::
-#
- #\begin{Bmatrix}\vec{\dot\sig}_u \\ \vec{\dot\sig}_p\end{Bmatrix}=\begin{bmatrix}\mat{K}_{uu} & \mat{K}_{up} \\ \mat{K}_{pu} & \mat{K}_{pp}\end{bmatrix}\begin{Bmatrix}\vec{\dot\eps}_u \\ \vec{\dot\eps}_p\end{Bmatrix}
-#
-#Since we can only apply strain (via the :yref:`velocity gradient<Cell.velGrad>`), we must evaluate
-#
-#.. math :: \vec{\dot\eps}_u=\mat{K}_{uu}^{-1}(\vec{\dot\sig}_u-\mat{K}_{up}\vec{\dot\eps}_p)
-#
-#The :yref:`velocity gradient <Cell.velGrad>` is then assembled from $\vec{\eps}_u$ and $\vec{\eps}_p$. If magnitude of any component of $\nabla\vec{v}$ exceeds :yref:`maxStrainRate<Peri3dController.maxStrainRate>`, the whole tensor is scaled accordingly.
-#
-#'''
+
+wrapper.Ig2_Sphere_Sphere_L3Geom_Inc.__doc__=r'''Functor for computing incrementally configuration of 2 :yref:`Spheres<Sphere>` stored in :yref:`L3Geom`; the configuration is positioned in global space by local origin $\vec{c}$ (contact point) and rotation matrix $\mat{T}$ (orthonormal transformation matrix), and its degrees of freedom are local displacement $\vec{u}$ (in one normal and two shear directions); with :yref:`Ig2_Sphere_Sphere_L6Geom_Inc` and :yref:`L6Geom`, there is additionally $\vec{\phi}$. The first row of $\mat{T}$, i.e. local $x$-axis, is the contact normal noted $\vec{n}$ for brevity. Additionally, quasi-constant values of $\vec{u}_0$ (and $\vec{\phi}_0$) are stored as shifted origins of $\vec{u}$ (and $\vec{\phi}$); therefore, current value of displacement is always $\curr{\vec{u}}-\vec{u}_0$.
+
+Suppose two spheres with radii $r_i$, positions $\vec{x}_i$, velocities $\vec{v}_i$, angular velocities $\vec{\omega}_i$.
+
+When there is not yet contact, it will be created if $u_N=|\curr{\vec{x}}_2-\curr{\vec{x}}_1|-|f_d|(r_1+r2)<0$, where $f_d$ is :yref:`distFactor<Ig2_Sphere_Sphere_L3Geom_Inc.distFactor>` (sometimes also called ``interaction radius''). If $f_d>0$, then $\vec{u}_{0x}$ will be initalized to $u_N$, otherwise to 0. In another words, contact will be created if spheres enlarged by $|f_d|$ touch, and the ``equilibrium distance'' (where $\vec{u}_x-\vec{u}-{0x}$ is zero) will be set to the current distance if $f_d$ is positive, and to the geometrically-touching distance if negative.
+
+Local axes (rows of $\mat{T}$) are initially defined as follows:
+
+* local $x$-axis is $\vec{n}=\vec{x}_l=\normalized{\vec{x}_2-\vec{x}_1}$;
+* local $y$-axis positioned arbitrarily, but in a deterministic manner: aligned with the $xz$ plane (if $\vec{n}_y<\vec{n}_z$) or $xy$ plane (otherwise);
+* local $z$-axis $\vec{z}_l=\vec{x}_l\times\vec{y}_l$.
+
+If there has already been contact between the two spheres, it is updated to keep track of rigid motion of the contact (one that does not change mutual configuration of spheres) and mutual configuration changes. Rigid motion transforms local coordinate system and can be decomposed in rigid translation (affecting $\vec{c}$), and rigid rotation (affecting $\mat{T}$), which can be split in rotation $\vec{o}_r$ perpendicular to the normal and rotation $\vec{o}_t$ (``twist'') parallel with the normal:
+
+.. math:: \pprev{\vec{o}_r}=\prev{\vec{n}}\times\curr{\vec{n}}.
+
+Since velocities are known at previous midstep ($t-\Dt/2$), we consider mid-step normal
+
+.. math:: \pprev{\vec{n}}=\frac{\prev{\vec{n}}+\curr{\vec{n}}}{2}.
+
+For the sake of numerical stability, $\pprev{\vec{n}}$ is re-normalized after being computed, unless prohibited by :yref:`approxMask<Ig2_Sphere_Sphere_L3Geom_Inc.approxMask>`. If :yref:`approxMask<Ig2_Sphere_Sphere_L3Geom_Inc.approxMask>` has the appropriate bit set, the mid-normal is not compute, and we simply use $\pprev{\vec{n}}\approx\prev{\vec{n}}$.
+
+Rigid rotation parallel with the normal is
+
+.. math:: \pprev{\vec{o}_t}=\pprev{\vec{n}}\left(\pprev{\vec{n}}\cdot\frac{\pprev{\vec{\omega}}_1+\pprev{\vec{\omega}}_2}{2}\right)\Dt.
+
+*Branch vectors* $\vec{b}_1$, $\vec{b}_2$ (connecting $\curr{\vec{x}}_1$, $\curr{\vec{x}}_2$ with $\curr{\vec{c}}$ are computed depending on :yref:`noRatch<Ig2_Sphere_Sphere_L3Geom_Inc.noRatch>` (see :yref:`here<Ig2_Sphere_Sphere_ScGeom.avoidGranularRatcheting>`).
+
+.. math::
+ :nowrap:
+
+ \begin{align*}
+ \vec{b}_1&=\begin{cases} r_1 \curr{\vec{n}} & \mbox{with \texttt{noRatch}} \\ \curr{\vec{c}}-\curr{\vec{x}}_1 & \mbox{otherwise} \end{cases} \\
+ \vec{b}_2&=\begin{cases} -r_2\curr{\vec{n}} & \mbox{with \texttt{noRatch}} \\ \curr{\vec{c}}-\curr{\vec{x}}_2 & \mbox{otherwise} \end{cases} \\
+ \end{align*}
+
+Relative velocity at $\curr{\vec{c}}$ can be computed as
+
+.. math:: \pprev{\vec{v}_r}=(\pprev{\vec{\tilde{v}}_2}+\vec{\omega}_2\times\vec{b}_2)-(\vec{v}_1+\vec{\omega}_1\times\vec{b}_1)
+
+where $\vec{\tilde{v}}_2$ is $\vec{v}_2$ without mean-field velocity gradient in periodic boundary conditions (see :yref:`Cell.homoDeform`). In the numerial implementation, the normal part of incident velocity is removed (since it is computed directly) with $\pprev{\vec{v}_{r2}}=\pprev{\vec{v}_r}-(\pprev{\vec{n}}\cdot\pprev{\vec{v}_r})\pprev{\vec{n}}$.
+
+Any vector $\vec{a}$ expressed in global coordinates transforms during one timestep as
+
+.. math:: \curr{\vec{a}}=\prev{\vec{a}}+\pprev{\vec{v}_r}\Dt-\prev{\vec{a}}\times\pprev{\vec{o}_r}-\prev{\vec{a}}\times{\pprev{\vec{t}_r}}
+
+where the increments have the meaning of relative shear, rigid rotation normal to $\vec{n}$ and rigid rotation parallel with $\vec{n}$. Local coordinate system orientation, rotation matrix $\mat{T}$, is updated by rows, i.e.
+
+.. math:: \curr{\mat{T}}=\begin{pmatrix} \curr{\vec{n}_x} & \curr{\vec{n}_y} & \curr{\vec{n}_z} \\ \multicolumn{3}{c}{\prev{\mat{T}_{1,\bullet}}-\prev{\mat{T}_{1,\bullet}}\times\pprev{\vec{o}_r}-\prev{\mat{T}_{1,\bullet}}\times\pprev{\vec{o}_t}} \\ \multicolumn{3}{c}{\prev{\mat{T}_{2,\bullet}}-\prev{\mat{T}_{2,\bullet}}\times\pprev{\vec{o}_r}-\prev{\mat{T}_{,\bullet}}\times\pprev{\vec{o}_t}} \end{pmatrix}
+
+This matrix is re-normalized (unless prevented by :yref:`approxMask<Ig2_Sphere_Sphere_L3Geom_Inc.approxMask>`) and mid-step transformation is computed using quaternion spherical interpolation as
+
+.. math:: \pprev{\mat{T}}=\mathrm{Slerp}\,\left(\prev{\mat{T}};\curr{\mat{T}};t=1/2\right).
+
+Depending on :yref:`approxMask<Ig2_Sphere_Sphere_L3Geom_Inc.approxMask>`, this computation can be avoided by approximating $\pprev{\mat{T}}=\prev{\mat{T}}$.
+
+Finally, current displacement is evaluated as
+
+.. math:: \curr{\vec{u}}=\prev{u}+\pprev{\mat{T}}\pprev{\vec{v}_r}\Dt.
+
+For the normal component, non-incremental evaluation is preferred, giving
+
+.. math:: \curr{\vec{u}_x}=|\curr{\vec{x}_2}-\curr{\vec{x}_1}|-(r_1+r_2)
+
+If this functor is called for :yref:`L6Geom`, local rotation is updated as
+
+.. math:: \curr{\vec{\phi}}=\prev{\vec{\phi}}+\pprev{\mat{T}}\Dt(\vec{\omega}_2-\vec{\omega}_1)
+
+.. note: TODO: ``distFactor`` is not yet implemented as described above; some formulas mix values at different times, should be checked carefully.
+
+'''
=== modified file 'py/tests/pbc.py'
--- py/tests/pbc.py 2010-12-10 20:20:16 +0000
+++ py/tests/pbc.py 2010-12-11 23:01:30 +0000
@@ -6,7 +6,7 @@
'''
import unittest
-import random
+import random,math
from yade.wrapper import *
from miniEigen import *
from yade._customConverters import *
@@ -32,20 +32,20 @@
O.cell.homoDeform=3
O.dt=0 # do not change positions with dt=0 in NewtonIntegrator, but still update velocities from velGrad
def testHomotheticResizeVel(self):
- "PBC: (homoDeform==3) homothetic cell deformation adjusts particle velocity"
+ "PBC: homothetic cell deformation adjusts particle velocity (homoDeform==3)"
O.dt=1e-5
O.step()
s1=O.bodies[1].state
self.assertAlmostEqual(s1.vel[2],self.initVel[2]+self.initPos[2]*O.cell.velGrad[2,2])
def testHomotheticResizePos(self):
- "PBC: (homoDeform==1) homothetic cell deformation adjusts particle position"
+ "PBC: homothetic cell deformation adjusts particle position (homoDeform==1)"
O.cell.homoDeform=1
O.step()
s1=O.bodies[1].state
self.assertAlmostEqual(s1.vel[2],self.initVel[2])
self.assertAlmostEqual(s1.pos[2],self.initPos[2]+self.initPos[2]*O.cell.velGrad[2,2]*O.dt)
def testScGeomIncidentVelocity(self):
- "PBC: (homoDeform==3) ScGeom computes incident velocity correctly"
+ "PBC: ScGeom computes incident velocity correctly (homoDeform==3)"
O.step()
O.engines=[InteractionLoop([Ig2_Sphere_Sphere_ScGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[])]
i=utils.createInteraction(0,1)
@@ -53,7 +53,7 @@
self.assertEqual(self.initVel,i.geom.incidentVel(i,avoidGranularRatcheting=False))
self.assertAlmostEqual(self.relDist[1],1-i.geom.penetrationDepth)
def testScGeomIncidentVelocity_homoPos(self):
- "PBC: (homoDeform==1) ScGeom computes incident velocity correctly"
+ "PBC: ScGeom computes incident velocity correctly (homoDeform==1)"
O.cell.homoDeform=1
O.step()
O.engines=[InteractionLoop([Ig2_Sphere_Sphere_ScGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[])]
@@ -61,8 +61,23 @@
self.assertEqual(self.initVel,i.geom.incidentVel(i,avoidGranularRatcheting=True))
self.assertEqual(self.initVel,i.geom.incidentVel(i,avoidGranularRatcheting=False))
self.assertAlmostEqual(self.relDist[1],1-i.geom.penetrationDepth)
+ def testL3GeomIncidentVelocity(self):
+ "PBC: L3Geom computes incident velocity correctly (homoDeform==3)"
+ O.step()
+ O.engines=[ForceResetter(),InteractionLoop([Ig2_Sphere_Sphere_L3Geom_Inc()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_L3Geom_FrictPhys_ElPerfPl(noBreak=True)]),NewtonIntegrator()]
+ i=utils.createInteraction(0,1)
+ O.dt=1e-10; O.step() # tiny timestep, to not move the normal too much
+ self.assertAlmostEqual(self.initVel.norm(),(i.geom.u/O.dt).norm())
+ def testL3GeomIncidentVelocity_homoPos(self):
+ "PBC: L3Geom computes incident velocity correctly (homoDeform==1)"
+ O.cell.homoDeform=1; O.step()
+ O.engines=[ForceResetter(),InteractionLoop([Ig2_Sphere_Sphere_L3Geom_Inc()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_L3Geom_FrictPhys_ElPerfPl(noBreak=True)]),NewtonIntegrator()]
+ i=utils.createInteraction(0,1)
+ O.dt=1e-10; O.step()
+ self.assertAlmostEqual(self.initVel.norm(),(i.geom.u/O.dt).norm())
+ #self.assertAlmostEqual(self.relDist[1],1-i.geom.penetrationDepth)
def testKineticEnergy(self):
- "PBC: (homoDeform==3) utils.kineticEnergy considers only fluctuation velocity, not the velocity gradient"
+ "PBC: utils.kineticEnergy considers only fluctuation velocity, not the velocity gradient (homoDeform==3)"
O.step() # updates velocity with homotheticCellResize
# ½(mv²+ωIω)
# #0 is still, no need to add it; #1 has zero angular velocity
@@ -70,7 +85,7 @@
Ek=.5*O.bodies[1].state.mass*self.initVel.squaredNorm()
self.assertAlmostEqual(Ek,utils.kineticEnergy())
def testKineticEnergy_homoPos(self):
- "PBC: (homoDeform==1) utils.kineticEnergy considers only fluctuation velocity, not the velocity gradient"
+ "PBC: utils.kineticEnergy considers only fluctuation velocity, not the velocity gradient (homoDeform==1)"
O.cell.homoDeform=1; O.step()
self.assertAlmostEqual(.5*O.bodies[1].state.mass*self.initVel.squaredNorm(),utils.kineticEnergy())
