yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2274
committer: Václav Šmilauer <eudoxos@xxxxxxxx>
branch nick: trunk
timestamp: Thu 2010-06-03 23:01:34 +0200
message:
  1. Update docs for Peri3dController
  2. Fix compiler warning
modified:
  doc/sphinx/conf.py
  gui/qt3/SimulationController.cpp
  pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.hpp
  py/_extraDocs.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 'doc/sphinx/conf.py'
--- doc/sphinx/conf.py	2010-05-24 21:14:57 +0000
+++ doc/sphinx/conf.py	2010-06-03 21:01:34 +0000
@@ -373,10 +373,31 @@
 \usepackage{euler} % must be loaded before fontspec for the whole doc (below); this must be kept for pngmath, however
 \usepackage{amsmath}
 \usepackage{amsbsy}
-\let\mat\boldsymbol
-\global\let\vec\boldsymbol
+
+% symbols
+\let\mat\boldsymbol % matrix
+\let\vec\boldsymbol % vector
+\let\tens\boldsymbol % tensor
+
+% algorithm complexity
+\def\bigO#1{\ensuremath{\mathcal{O}(#1)}}
+
+% variants for greek symbols
+\let\epsilon\varepsilon
 \let\theta\vartheta
 \let\phi\varphi
+
+% shorthands
+\let\sig\sigma
+\let\eps\epsilon
+
+% variables at different points of time 
+\def\prev#1{#1^-}
+\def\pprev#1{#1^\ominus}
+\def\curr#1{#1^{\circ}}
+\def\nnext#1{#1^\oplus}
+\def\next#1{#1^+}
+
 '''
 
 pngmath_latex_preamble=r'\usepackage[active]{preview}\usepackage{amsmath}\usepackage{amssymb}'+my_latex_preamble

=== modified file 'gui/qt3/SimulationController.cpp'
--- gui/qt3/SimulationController.cpp	2010-05-30 13:41:00 +0000
+++ gui/qt3/SimulationController.cpp	2010-06-03 21:01:34 +0000
@@ -394,7 +394,6 @@
 		unsigned int hours	= duration.hours();
 		unsigned int minutes 	= duration.minutes();
 		unsigned int seconds	= duration.seconds();
-		unsigned int mseconds	= duration.fractional_seconds()/1000;
 		unsigned int days 	= hours/24;
 		hours			= hours-24*days;
 

=== modified file 'pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.hpp'
--- pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.hpp	2010-05-26 17:49:15 +0000
+++ pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.hpp	2010-06-03 21:01:34 +0000
@@ -83,7 +83,7 @@
 
 		virtual void action();
 		void update();
-	YADE_CLASS_BASE_DOC_ATTRS(Peri3dController,BoundaryController,"Experimental controller of full strain/stress tensors on periodic cell. Stress and strain tensors are computed using formulas derived in [Kuhl2001]_, in particular equations (33) and (35).",
+	YADE_CLASS_BASE_DOC_ATTRS(Peri3dController,BoundaryController,"Experimental controller of full strain/stress tensors on periodic cell. Detailed documentation is in py/_extraDocs.py.",
 		((Matrix3r,strain,Matrix3r::Zero(),"Current deformation tensor |yupdate|"))
 		((Matrix3r,stress,Matrix3r::Zero(),"Current stress tensor |yupdate|"))
 		((Matrix3r,goal,Matrix3r::Zero(),"Goal state."))

=== modified file 'py/_extraDocs.py'
--- py/_extraDocs.py	2010-06-03 07:31:33 +0000
+++ py/_extraDocs.py	2010-06-03 21:01:34 +0000
@@ -30,26 +30,77 @@
 
 
 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.
+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::
-
-	\begin{align}
-		\tens{K}&=\frac{1}{V}\sum_c {d_0^c}^2\left[K_N^c \vec{N}^c\otimes\vec{N}^c+K_T^c{\vec{T}^c}^T\cdot\vec{T}^c
-
-
-
-
-
-
-
-
-
-
-
+	: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 should be obtained via
+
+.. math:: \mat{\eps}=\mat{U}\log \mat{P}
+
+but as matrix logarithm is not yet implemented in Eigen, we approximate this (which works only for "small strain" scenarios) as
+
+.. math:: \mat{\eps}=\mat{U}(\mat{P}-\mat{I}_3)
+
+where $\mat{I}_3$ is 3×3 identity matrix.
+
+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}$ (dealing with large strains, we should use the rate formulation; supposing $\mat{K}$ is quasi-constant in time (linearization), we have $\vec{\dot\sig}=\mat{K}\vec{\dot\eps}$, where all rate terms differ from their counterparts by the same factor $\Delta t$. For that reason, it is only at the end that the rate form is assumed in the proper sense) is decomposed in 2 parts:
+
+* $p$ denotes the part where strain component $\vec{\eps}_i$ is *prescribed*;
+* $u$ denotes the part where strain component $\vec{\eps}_i$ is *unprescribed* (and $\vec{\sig}_i$ is prescribed).
+
+The equation then is re-arranged as
+
+.. math:: 
+
+	\begin{Bmatrix}\vec{\sig}_u \\ \vec{\sig}_p\end{Bmatrix}=\begin{bmatrix}\mat{K}_{uu} & \mat{K}_{up} \\ \mat{K}_{pu} & \mat{K}_{pp}\end{bmatrix}\begin{Bmatrix}\vec{\eps}_u \\ \vec{\eps}_p\end{Bmatrix}
+
+Since we can only apply strain (via the :yref:`velocity gradient<Cell.velGrad>`), we must evaluate 
+
+.. math :: \vec{\eps}_u=\mat{K}_{uu}^{-1}(\vec{\sig}_u-\mat{K}_{up}\vec{\eps}_p)
+
+The goal strain tensor $\tens{\eps}^g$ is the assembled from $\vec{\eps}_u$ and $\vec{\eps}_p$ and is used to prescribe the current :yref:`velocity gradient <Cell.velGrad>`
+
+.. math:: \nabla\vec{v}=\frac{\tens{\eps}^g-\tens{\eps}}{\Delta t}.
+
+If magnitude of any component of $\nabla\vec{v}$ exceeds :yref:`maxStrainRate<Peri3dController.maxStrainRate>`, the whole tensor is scaled accordingly.
 
 '''