Code::Blocks + GLUT + FLTK

Code::Blocks provides project templates for both GLUT and FLTK, but no for GLUT + FLTK. Here is a procedure to build a simple GLUT + FLTK program in Code::Blocks.


1.  Build up a new project for GLUT by using the Code::Blocks' template.
2.  Add search directories setting to the include and library directories of FLTK.
3.  Modify #include <GL/glut.h> in main.cpp to #include<FL/glut.H>
4.  Modify the Linker settings as follows
   1)  Add other linker options:  -mwindows
   2)  Link libraries as: fltk_images, fltk_gl, fltk, fltk_png, fltk_z, fltk_jpeg, glu32, opengl32, ole32, uuid, comctl32. Pay attention to that the sequential order of those libraries is important.


The project would be rebuilt successfully. You can then add FLTK functions as you like.

How to install GLUT for mingw

There are some binary distribution of glut in windows available, such as here. But it is not sure if they are available for your system, the version of your gcc compiler on your mingw32 or mingw64 etc. On the other hand, it is not easy to compile its sourcecode in mingw. However, freeglut provides cmake support for cross platform compiling. It is therefore suggested that freeglut, which provides binary distribution also, be adopted.


1. Compile freeglut by using cmake.


2. Keeping the include, library  and dll file of freeglut available and compile like:
    g++ -Os -Wall test.cc -o test -mwindows -lfreeglut -lopengl32 -lglu32


That's all.

如何将windows的dll库文件转换为mingw库文件

1  下载gendef
http://sourceforge.net/projects/mingw/files/MinGW/Extension/gendef/gendef-1.0.1346/


2.  执行gendef dll库文件名产生def文件
   如  gendef msmpi.dll 将产生 msmpi.def


3   执行dlltool -d msmpi.def -l libmpi.def -l libmsmpi.a -D msmpi.dll 将库文件名产生libmsmpi.a


如果要在64位机上使用msmpi，除此以外还需
1)  在/Inc/mpi.h文件的第128行，加上
     #include <stdint.h>
2） 如果使用非intel编译器，如gfortran，还需将/Inc/mpif.h文件用用mpich中的相关文件替代。

Build Trilinos in VS+intel compiler

Following software are used:
1). Trilinos 11.14.1
2). Visual studio 2012
3). Intel compiler Composer XE2013
4). cmake 3.2.1


1.  Begin with intel Command Prompt Visual Studio 2012 mode
2.  Type in : cmake-gui to open cmake gui
3.  First configure in cmake
     Specif the generate for this project: " Visual studio 11 2012 Win64"
     Speify native compiler: Here you set CMAKE_C_COMPILER.  CAMKE_CXX_COMPILER and CMAKE_FORTRAN_COMPILER to icl, ifort respectively.


4.  Afetr first configuration, Turn off Tilinos_ENABLE_OPTIONAL_PACKAGES because most libraries in trilinos cannot be compiled in windows. Turn on Trilinos_ENABLE_Epetra, Trilinos_ENABLE_Epetraext, Trilinos_ENABLE_ML. Turn on Triluinos_ENABLE_Fortran, TPL_ENABLE_MPI if needed.


5. Do configure


6. Libraries like BLAS, LAPACK, MPI maybe missing. Set
   BLAS_LIBRARIES_DIRS = /pahto/mkl/lib/intel64
   BLAS_LIBRARIE_MAKE = mkl_intel_lp64


   LAPACK_LIBRARIES_DIRS = /pahto/mkl/lib/intel64
   LAPACK_LIBRARIE_MAKE = mkl_intel_lp64


  MPI_BASE_DIR = C:/Program Files/Microsoft HPC Pack2012/
  ........
by hand.


7.  Do configure
8. Do generate. You would get MSVC solution files.
9. Double click Trillinos.sln.
10. Add a preprocessor definition ICL and compile.

Build Opensees in Visual studio 2012

SVN version of Opensees: 5841

OS:  Windows 7

Visual studio 2012

The compileration use the attached Microsoft Visual Studio Solution OpenseesVC6.sln directly.

1. Download and install ActiveTcl 8.5 for x86. ActiveTcl 8.6 deletes some class used by Opensees. 64bit ActiveTcl is also not available, at least without much modification in OpenseesVC6.sln.

2. You need add the position of Tcl include and library folder by hand for compilering some of the Opensees library.

3.  You need close SAFESEH in linker option for compilering some of the Opensees library.

4.  Following files are lost from the compiler list, you must add them by hand.

    Pileto3D.cpp,  FRPConfinedConcrete.cpp, YamamotoBiaxialHDR.cpp

5.  The file /win32/VC2005errno.cpp is also needed to be added in compile list otherwise you would get unrresolved _errono error.

That's all.

逆变基底不能随体变化吗？

有网友指出*两本教科书（参考文献[1]，[2]）上指出逆变基底不能随体变化。其原文如下

如果这一结论成立，一个很麻烦的结果是你将得不到逆变基矢量的随体微分（convected derivative）,而没有这一客观性微分整个涉及逆变基矢量的运算法则都要重写。另一个诡异的结果出自参考文献[1]，对应于上面的文字，作者认为参考构型的逆变基底是时间的函数！（图1中的$G^1(t), G^2(t)$: Has to be reformed every moment with a deformation.）

图1

那么问题出在哪里呢？

参考文献[1]，[2]的作者都一下面的图来解释逆变基矢量的几何意义, 即将协，逆变基矢都理解为印刻在变形体有长度方向的线段(图2)。当物体变形是，协变基矢$a_1,a_2$随变形体一起变动, 而逆变基矢$a^1,a^2$不能随变形体一起变动。两位作者都是这样解释的。

图2

  问题就出在上述逆变基矢量的几何解释上。对于数学概念，It often seems like there are two types of students of mathematics: those who prefer to learn by studying equations and following derivations, and those who prefer pictures（[10]). 学数学的往往是第一类，更由于几何表示逆变基底或协变矢量比较复杂，你很少能看到在数学书上看到逆变基底或协变矢量的几何表示, 这对希望得到几何图像的第二类型的人造成了不少困惑（如[3-6]）**。 学工程的则是第二类，你可以在非数学类的教科书中看到图二的解释。这样的解释也许对理解逆变基矢量有所帮助，但是要了解这一解释并不严格，只能用于帮助你理解而并不是严格的定义，否则会如参考文献[1]，[2]得出错误的结果。
  下面图3关于逆变矢量的几何解释出自文献[6-9]。在这里协变矢量被表示为一组面。 这一解释基于逆变矢量和协变矢量的对偶关系:一个逆（协）变矢量作用到协（逆）变矢量产生一个实数。这种几何解释也许不太好理解，但这里要强调的是逆变基底不能理解为上述作者所述附在变形体上的线段！在参考文献[6]中，逆变基底被解释为坐标超平面的法线，其随体变化则表征的该超平面的变化，是存在的。实际上，逆变矢量的convected derivative可以从其与协变矢量的对偶关系中推导出来。
  

图3 逆变矢量的几何解释

  如上所述，正是由于参考文献[1]，[2]的作者对逆变矢量几何意义的错误解释导致了逆变基底不随体变化的错误结论。其实我们并不需要对所有的数学概念画图，不信你画一个六维的矢量，四阶的张量试试😀

  最后介绍一个可以直观地表示逆变，协变矢量的软件GeoGeobra。也算是给喜欢画图的同学一个交代！

1.  Koichi Hashiguchi, Yuki Yamakawa:  Introduction to Finite Strain Theory for Continuum Elasto-Plasticity, Wiley, 2012
2.  黄克智:  张量分析, 清华大学出版社; 第2版 (2003年7月1日)
3.  Geometric imagine of one forms: http://mathoverflow.net/questions/26939/geometric-imagination-of-differential-forms
4. Physical and geometric interpretation of Differential form:  http://physics.stackexchange.com/questions/55751/physical-and-geometrical-interpretation-of-differential-forms
5. Help with geometric interpretation of 1-form: https://www.physicsforums.com/threads/help-with-geometric-interpretation-of-1-form.376475/
6. Are there pictorial examples that distinguish covariant and contravariant vectors?
https://arxiv.org/pdf/gr-qc/9807044.pdf
7  https://en.wikipedia.org/wiki/Covariance_and_contravariance_of_vectors
8  http://en.wikipedia.org/wiki/Linear_form
9.  B.F. Shutz: A First Course in General Relativity, Cambridge, 2009
10  D. Bachman: A Geometric Approach to Differential Form, Springer Science+Business Media, LLC 2006, 2012
11 https://www.geogebra.org/

注:
* 本文问题的提出出自http://forum.simwe.com/thread-1115839-1-1.html
** 张量分析中的逆变矢量，协变矢量也有称为矢量（vector）和余矢量或对偶矢量（covector，dual vector）的，在量子力学（狄拉克记号）中被称为ket和bra，在微分几何中则被称为切矢（tangent vector）和余切矢（cotagent），或矢量和1形式。

有限元接触解析算法（之二）

罚函数法，Lagrangian乘数法，递增Lagrangian乘数法， 内点法

3  渐增Lagrangian乘数法/内点法的实装方法

前述第1节讲述的四种方法中，罚函数法，Lagrangian乘数法的解可以从所得方程直接解出。而递增Lagrangian乘数法和内点法需要使用迭代法求解。

在本节中，将以一自由度的接触问题为例，又易到难地讲解采用递增Lagrangian乘数法和内点法的接触算法。由于接触问题的解析一般用增分法求解，而在第一节中我们给出的是全量公式。在这一节中给出的算法是基于增分法的，因此第一节中给出的公式也要相应改写。

为显示算法实现，从本节开始将添加一些用python语言编写的程序附件。在之后的算法实现中。将大量使用迭代计算。为了明示迭代步数，本文将以上下标n，k表示各迭代步。

3.1  考虑接触非线性的基于递增Lagrangian乘数法的实装方法

图1 一自由度的接触问题

S0:  初始状态的求解

  1） 初期接触状态的判定:无接触，由（0.4）式得到物体m的位移量$x_0=F/K$

  2)  接触状态的判定: 如$g=g_0+x_0>=0$计算结束。否则判定进入接触，继续下一步计算。

  以下进入递增Lagrangian乘数法求解。在这里我们令罚函数$\varepsilon_N=10000$,收束判定值TOL=1.e-10

在此，对应于增分解析，我们将(3.1)改写为
  $J=1/2K(x+x_0)^2 -F(x+x_0) +\lambda (x+g)+1/2\varepsilon_N (x+g)^2=0$
在这里x表示下一步的位移。相应地，式(3.4)改变为
  $ x=(F-\varepsilon_Ng_0-\lambda-Kx_0)/(K+\varepsilon_N)  $

S1: k=0, $\lambda_0=0$, $x_0=0$
S2: if $g+x_k<=TOL$
         goto S3
      else
         $ x_{k+1}=(F-\varepsilon_Ng_0-\lambda-Kx_0)/(K+\varepsilon_N)  $
         $\lambda_{k+1}=\lambda_{k}+ \varepsilon_N *(g+x_{k+1})$
         k=k+1;  goto S2
S3: end

计算源码从此下载。

计算过程如下。
iter,    disp  incr   ,        gap     ,       Lagrangian
  0,  0.000000000, -1.0000000e-01,  0.000000000
  1,  0.090909091, -9.0909091e-03, -9.090909091
  2,  0.099173554, -8.2644628e-04, -9.917355372
  3,  0.099924869, -7.5131480e-05, -9.992486852
  4,  0.099993170, -6.8301346e-06, -9.999316987
  5,  0.099999379, -6.2092132e-07, -9.999937908
  6,  0.099999944, -5.6447393e-08, -9.999994355
  7,  0.099999995, -5.1315812e-09, -9.999999487
  8,  0.100000000, -4.6650737e-10, -9.999999953
  9,  0.100000000, -4.2409770e-11, -9.999999996
Displacement= -0.100000000042
Final gap= -4.24097701401e-11
Contact force= -9.99999999576

3.2  考虑接触非线性的的基于内点法的实装方法

这里的算法在（1.4）节给出的primal-dual算法。计算源码从此下载。

计算过程如下。

iter, displacement, Lagrangian, barrier parameter

  1, -0.038196601, 16.180339887, 1.0000000e+00

  2, -0.075838015, 12.416198487, 3.0000000e-01

  3, -0.091690481, 10.830951895, 9.0000000e-02

  4, -0.097369211, 10.263078947, 2.7000000e-02

  5, -0.099196457, 10.080354318, 8.1000000e-03

  6, -0.099757588, 10.024241236, 2.4300000e-03

  7, -0.099927153, 10.007284693, 7.2900000e-04

  8, -0.099978135, 10.002186522, 2.1870000e-04

  9, -0.099993439, 10.000656057, 6.5610000e-05

 10, -0.099998032, 10.000196826, 1.9683000e-05

 11, -0.099999410, 10.000059049, 5.9049000e-06

 12, -0.099999823, 10.000017715, 1.7714700e-06

 13, -0.099999947, 10.000005314, 5.3144100e-07

 14, -0.099999984, 10.000001594, 1.5943230e-07

 15, -0.099999995, 10.000000478, 4.7829690e-08

 16, -0.099999999, 10.000000143, 1.4348907e-08

 17, -0.100000000, 10.000000043, 4.3046721e-09

 18, -0.100000000, 10.000000013, 1.2914016e-09

 19, -0.100000000, 10.000000004, 3.8742049e-10

 20, -0.100000000, 10.000000001, 1.1622615e-10

Displacement= -0.0999999999884

Final gap= 1.16226195335e-11

Contact force= 10.0000000012

3.3  考虑材料非线性和接触非线性的基于递增Lagrangian乘数法的实装方法

在前述算例中我们只考虑了接触非线性，而在实际应用中，接触非线性往往和材料非线性，几何非线性等并存。这一节将考虑非线性弹簧，其弹性系数K=100-1/30x。以此演示和其他非线性现象共存时接触解析的基于递增Lagrangian乘数法的实装方法。

此时的泛函J的变分为
$J=1/2K_0^2-1/30 x^3+ (\lambda+1/2\varepsilon_Nx) x -F x$
$\delta J=(K_0-0.1x) x\delta x+ (\lambda+\varepsilon_Nx)\delta x -(F-\varepsilon_Ng_0) \delta x=0$
基于变分$\delta x$的任意性，可得
$(K-0.1x+\varepsilon_N)x+\lambda-(F+\varepsilon_Ng_0)=0$
这是在这一节要解的方程。采用Newton迭代法，上式的增分形式为
$\Delta x^{n+1}=(F-(K_0-0.1x^n)x^n-\lambda-\varepsilon_N(x^n+g_0))/(K_0-0.2x^n+\varepsilon_N)$

S0:  初始状态的求解

  1） 初期接触状态的判定:无接触，由$(K+0.1x)x+\lambda=F$求解物体m的初始位移量$x_0$。由于该式非线性，我们需要使用迭代计算。

  2)  接触状态的判定: 如$g=g_0+x_0>=0$计算结束。否则判定进入接触，继续下一步计算。

  以下进入递增Lagrangian乘数法求解。在此我们使用两重迭代(nested augmented lagrangian algorithm)，分别处理材料非线性和接触非线性
S1: k=0, $\lambda_0=0$, $x_0^0=0$
S2: 外部循环，更新计算Lagrangian乘数
      if $g+x_k<=TOL$
         goto S5
      else
         goto S3
S3   内部循环，Lagrangian乘数不变，更新计算位移增量
       n = 0
       $df =F-\varepsilon_Ng_0-\lambda-(K_0-0.1 x_k^n) x_k^n$
       if   df<TOL
          goto S4
       else
          $ \Delta x_{k+1}^{n+1}=df/(K-0.2x_k^n+\varepsilon_N)  $
           $x_{k+1}^{n+1}= x_{k+1}^n + \Delta x_{k+1}^{n+1}$
          n = n+1
          goto S3
S4   $\lambda_{k+1}=\lambda_{k}+ \varepsilon_N *(g+x_{k+1})$
      k=k+1;  goto S2
S5   end

计算源码从此下载。

计算过程如下。

Initial displacement and contact gap -0.196152422706 0.1
iter, disp incr     ,      gap     , Lagrangian
  0, -0.196152423, -9.6152423e-02,  0.000000000
         internal loop, disp= 0.0870619064302
         internal loop, disp= 6.89056837473e-07
  1, -0.109089827, -9.0898272e-03, -9.089827219
         internal loop, disp= 0.0082633153912
         internal loop, disp= 6.20737540768e-09
  2, -0.100826506, -8.2650562e-04, -9.916332840
         internal loop, disp= 0.000751354971978
         internal loop, disp= 5.13203608441e-11
  3, -0.100075151, -7.5150597e-05, -9.991483437
         internal loop, disp= 6.83174816365e-05
         internal loop, disp= 4.24294993083e-13
  4, -0.100006833, -6.8331151e-06, -9.998316552
         internal loop, disp= 6.21180988284e-06
  5, -0.100000621, -6.2130523e-07, -9.998937857
         internal loop, disp= 5.64812673181e-07
  6, -0.100000056, -5.6492560e-08, -9.998994350
         internal loop, disp= 5.13559388072e-08
  7, -0.100000005, -5.1366210e-09, -9.998999486
         internal loop, disp= 4.66957053199e-09
  8, -0.100000000, -4.6705044e-10, -9.998999953
         internal loop, disp= 4.24583585173e-10
  9, -0.100000000, -4.2466849e-11, -9.998999996
Displacement= -0.100000000042
Final gap= -4.24668494814e-11
Contact force= -9.99899999575

3.4  摩擦接触问题的求解方法

     

图2 一自由度的摩擦接触问题

下面将以图2的问题为例，讨论解析摩擦接触约束的方法。与上节相同，本节也采用弹性系数K=100+1/10x的非线性弹簧。

在这里我们采用库仑摩擦定理
(4.1)$\begin{eqnarray}
f<=-\mu*mg*sign \dot x \\
x=0;\    if\   f<\mu*mg
\end{eqnarray}$
在$f<\mu*mg$是为粘结摩擦状态，否则称为滑动摩擦状态。

此时的泛函J为*1
(4.2)$\begin{eqnarray}
在粘结摩擦状态下:  J=1/2K_0 x^2 +0.1/3 x^3 -F x +\lambda x+1/2\varepsilon_T x^2 \\
在滑动摩擦状态下:  J=1/2K_0 x^2 +0.1/3 x^3 -F x +\mu*mg x
\end{eqnarray}$
在这里，我们使用递增Lagrangian乘数法处理条件$x=0;\    if\   f<\mu*mg$，该泛函与式（3.1）相同。滑动摩擦状态下的泛函则对应于式（0.1）。

由于要考虑粘结和滑动的状态变化，其解法要相应地复杂一些。这一问题的解法如下:

S1: k=0, $\lambda_0=0$, $x_0^0=0$
S2: 外部循环，更新计算Lagrangian乘数
     if $\Delta x<=TOL$
         goto S5
      else
         goto S3
S3  内部循环，Lagrangian乘数不变，更新计算位移增量
      n=0
      $fric=\lambda_k^n$
       if 粘结摩擦
          $df = F-\lambda_k-\varepsilon_T x -(K_0+0.1x_k)x_k$
          $\Delta x_{n+1} =df/(K_0+0.2x_n+\varepsilon_T)$
          $fric = fric+\varepsilon_T x$
          if $fric>\mu*mg$
              $fric = \mu*mg$进入滑动摩擦
       else
          $df = F-fric-(K_0+0.1x_k)x_k$
          $\Delta x_{n+1} =df/(K_0+0.2x_n)$
          $fric = \mu*mg$
       $ x=x+\Delta x$
        if df<TOL
           goto S4
       else
           n=n+1;  goto S3
S4   if 粘结摩擦
          $\lambda_{k+1}=\lambda_{k}+ \varepsilon_T *x_{k+1}$
      else
           $\lambda_{k+1}=\mu*mg$
      k=k+1;  goto S2
S3   end

在这里我们显示的是一维问题。在2维，3维计算中，我们需要同时考虑3.3节所示不可贯入条件和本节所示摩擦条件，此时我们需要结合3.3和3.4节的内容。其详细将在后续博文中论及。

为考虑粘结摩擦到滑动摩擦变化的处理方法，所付计算源码（从此下载）考虑了一个3步加载的问题，其计算结果如下

-----Step: 1
External Force= 10.0 0.0
Iter, Displacement , Lagrangian
         internal loop, disp= 0.00909090909091 0.00909090909091
         internal loop, disp= 0.00909090157777 0.00909090157777
         internal loop, disp= 0.00909090157777 0.00909090157777
  1,  0.009090902,  9.090901578
         internal loop, disp= 0.000826459258203 0.000826459258203
         internal loop, disp= 0.000826453049022 0.000826453049022
         internal loop, disp= 0.000826453049022 0.000826453049022
  2,  0.000826453,  9.917354627
         internal loop, disp= 7.51322082622e-05 7.51322082622e-05
         internal loop, disp= 7.51321569456e-05 7.51321569456e-05
         internal loop, disp= 7.51321569456e-05 7.51321569456e-05
  3,  0.000075132,  9.992486784
         internal loop, disp= 6.830197019e-06 6.830197019e-06
         internal loop, disp= 6.83019659489e-06 6.83019659489e-06
         internal loop, disp= 6.83019659489e-06 6.83019659489e-06
  4,  0.000006830,  9.999316980
         internal loop, disp= 6.20926970885e-07 6.20926970885e-07
         internal loop, disp= 6.20926967378e-07 6.20926967378e-07
  5,  0.000000621,  9.999937907
         internal loop, disp= 5.64479061917e-08 5.64479061917e-08
         internal loop, disp= 5.64479061623e-08 5.64479061623e-08
  6,  0.000000056,  9.999994355
         internal loop, disp= 5.13162783284e-09 5.13162783284e-09
         internal loop, disp= 5.13162783284e-09 5.13162783284e-09
  7,  0.000000005,  9.999999487
         internal loop, disp= 4.66511621933e-10 4.66511621933e-10
         internal loop, disp= 4.66511620318e-10 4.66511620318e-10
  8,  0.000000000,  9.999999953
         internal loop, disp= 4.24101488013e-11 4.24101488013e-11
         internal loop, disp= 4.24101488013e-11 4.24101488013e-11
  9,  0.000000000,  9.999999996
         internal loop, disp= 3.85546854484e-12 3.85546854484e-12
         internal loop, disp= 3.85546692997e-12 3.85546692997e-12
 10,  0.000000000, 10.000000000
Total Displacement & Displacement of current step= 3.85546692997e-12 3.85546692997e-12
frictional force= 9.99999999961
Spring force= 3.85546692997e-10
Reaction force= 10.0

-----Step: 2
External Force= 20.0 9.99999999961
Iter, Displacement , Lagrangian
   From stick to slip state: 19.0909090905 > 10.0
         internal loop, disp= 0.00909090909476 0.00909090909091
         internal loop, disp= 0.0999982644944 0.0999982644905
         internal loop, disp= 0.0999900019996 0.0999900019957
         internal loop, disp= 0.0999900019995 0.0999900019956
 11,  0.099990002, 10.000000000
         internal loop, disp= 0.0999900019995 0.0999900019956
 12,  0.099990002, 10.000000000
Total Displacement & Displacement of current step= 0.0999900019995 0.0999900019956
frictional force= 10.0
Spring force= 10.0
Reaction force= 20.0

-----Step: 3
External Force= 10.0 10.0
Iter, Displacement , Lagrangian
   From slip to stick state: -89.980005998 < 10.0
         internal loop, disp= 9.99600149944e-06 -0.099980005998
         internal loop, disp= 0.0909000016525 -0.00909000034697
         internal loop, disp= 0.0908992506656 -0.00909075133393
         internal loop, disp= 0.0908992506656 -0.00909075133393
 13,  0.090899251, -9.090751334
         internal loop, disp= 0.0991634334768 -0.000826568522726
         internal loop, disp= 0.0991634272681 -0.000826574731406
         internal loop, disp= 0.0991634272681 -0.000826574731406
 14,  0.099163427, -9.917326065
         internal loop, disp= 0.0999148452942 -7.51567052527e-05
         internal loop, disp= 0.0999148452429 -7.51567565817e-05
         internal loop, disp= 0.0999148452429 -7.51567565817e-05
 15,  0.099914845, -9.992482822
         internal loop, disp= 0.0999831683259 -6.83367359656e-06
         internal loop, disp= 0.0999831683255 -6.83367402093e-06
         internal loop, disp= 0.0999831683255 -6.83367402093e-06
 16,  0.099983168, -9.999316496
         internal loop, disp= 0.0999893806435 -6.21356025034e-07
         internal loop, disp= 0.0999893806435 -6.21356028544e-07
 17,  0.099989381, -9.999937852
         internal loop, disp= 0.0999899455023 -5.64971807554e-08
         internal loop, disp= 0.0999899455023 -5.64971807833e-08
 18,  0.099989946, -9.999994349
         internal loop, disp= 0.0999899968625 -5.13704107356e-09
         internal loop, disp= 0.0999899968625 -5.1370410737e-09
 19,  0.099989997, -9.999999486
         internal loop, disp= 0.0999900015324 -4.67088634118e-10
         internal loop, disp= 0.0999900015324 -4.67088632216e-10
 20,  0.099990002, -9.999999953
         internal loop, disp= 0.099990001957 -4.24703217975e-11
         internal loop, disp= 0.099990001957 -4.24703219933e-11
 21,  0.099990002, -9.999999996
         internal loop, disp= 0.0999900019956 -3.86164012211e-12
         internal loop, disp= 0.0999900019956 -3.8616413703e-12
 22,  0.099990002, -10.000000000
Total Displacement & Displacement of current step= 0.0999900019956 -3.8616413703e-12
frictional force= 3.86238596661e-10
Spring force= 9.99999999961
Reaction force= 10.0

第一加载步处于粘结摩擦状态，可以看出此时的位移为零，但出现摩擦力。第二加载步则进入滑动摩擦，可以看出明显的位移，同时摩擦力满足库仑摩擦定理。在第3加载步，物体滑动方向逆转，再次进入粘结摩擦状态。在这一步，位移增分为零，摩擦力回零。

下章预告: 动力学数值解析中接触条件的处理方法

*1 对于摩擦状态变化的问题，没有统一的泛函表达式，需要使用虚位移原理公式。这一内容将在后续博文中论及。