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.
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 逆变矢量的几何解释
最后介绍一个可以直观地表示逆变，协变矢量的软件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表示各迭代步。
在本节中，将以一自由度的接触问题为例，又易到难地讲解采用递增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) $
在此，对应于增分解析，我们将(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
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)$
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
此时的泛函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.1S1: 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
计算源码从此下载。
计算过程如下。
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 一自由度的摩擦接触问题
在这里我们采用库仑摩擦定理
(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节的内容。其详细将在后续博文中论及。
-----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 对于摩擦状态变化的问题，没有统一的泛函表达式，需要使用虚位移原理公式。这一内容将在后续博文中论及。
有限元接触解析方法（之一）
罚函数法，Lagrangian乘数法，递增Lagrangian乘数法*1， 内点法
0 本文的用意和写法
该系列博文面对对有限元解析有所了解，但未接触过接触解析的初学者。本文以实用性为第一目标，试图在概念准确的前提下追求直观描述，特别强调给出某些实用算法的细节，并附有相应程序。由于这些算法细节都经过本人实装测试，其有效性是有一定的保证的。希望读完此系列博文可以自己动手组装自己的解析程序，并希望能够指出其中问题或错误，共同提高。本文的解析对象是泛函的极值问题，这是本文的起点。对泛函的极值问题不了解的读者可以跳过泛函的变分运算步骤，着重于其后的物理概念和算法的讨论。
1 考虑约束条件的最优化算法
泛函的极值问题的求解从数值计算方法的角度来说是最优化算法，而接触问题为其约束条件。罚函数法，Lagrangian乘数法，递增Lagrangian乘数法， 内点法(Interior point method, barrier method或内部罚函数法)都是标准的考虑约束条件的最优化算法。这一节将以图1的问题为例解释这几种解法。
图1 一自由度的接触问题* 2
图1所示问题如没有接触约束，其泛函J为
(0.1) $J=1/2K x^2 -F x$
这里我们将m的位移记为x，外力F<0。该泛函的极值条件为其变分为零，即
(0.2) $\delta J=K x\delta x -F \delta x=0$
基于变分$\delta x$的任意性，要求
(0.3) $K x -F=0$
即得
(0.4) $ x=F/K $
在此我们假设F/K足够大使得$|x|>g_0$，但是接触条件要求物体m不能进入接触面, 即$g=g_0+x>=0$。如g<0, 可以用下述方法导入接触约束。
1.1 罚函数法
采用罚函数法考虑接触约束，其泛函J为(1.1) $J=1/2K x^2 -F x +1/2\varepsilon_N (x+g_0)^2$
这里$\varepsilon_N>0$为罚函数。该泛函的极值条件为其变分为零，即
(1.2) $\delta J=(K+\varepsilon_N) x\delta x -(F-\varepsilon_Ng_0) \delta x=0$
基于变分$\delta x$的任意性，要求
(1.3) $(K+\varepsilon_N) x -(F-\varepsilon_Ng_0)=0$
即得
(1.4) $ x=(F-\varepsilon_Ng_0)/(K+\varepsilon_N) $
当罚函数$\varepsilon_N$足够大时，从该公式可以得到逼近于精确解$x=-g_0$的近似解。
1.2 Lagrangian乘数法
采用Lagrangian乘数法考虑接触约束，其泛函J为
(2.1) $J=1/2K x^2 -F x +\lambda (x+g_0)$
这里$\lambda$为Lagrangian乘数,这是一个待定自变量。泛函J的极值条件为其变分为零，即
(2.2) $\delta J=(Kx+\lambda-F)\delta x +\delta \lambda (x+g_0)=0$
基于变分$\delta x$和$\delta \lambda$的任意性，要求
(2.3) $\begin{eqnarray}
Kx+\lambda=F \\
x=-g_0
\end{eqnarray}$
即得
(2.4) $ x=-g_0; \lambda=F+Kg_0 $
这是该问题的精确解。
(2.1) $J=1/2K x^2 -F x +\lambda (x+g_0)$
这里$\lambda$为Lagrangian乘数,这是一个待定自变量。泛函J的极值条件为其变分为零，即
(2.2) $\delta J=(Kx+\lambda-F)\delta x +\delta \lambda (x+g_0)=0$
基于变分$\delta x$和$\delta \lambda$的任意性，要求
(2.3) $\begin{eqnarray}
Kx+\lambda=F \\
x=-g_0
\end{eqnarray}$
即得
(2.4) $ x=-g_0; \lambda=F+Kg_0 $
这是该问题的精确解。
1.3 递增Lagrangian乘数法
采用递增Lagrangian乘数法考虑接触约束，其泛函J为*3
(3.1) $J=1/2K x^2 -F x +\lambda (x+g_0)+1/2\varepsilon_N (x+g_0)^2$
注意到该泛函既包括Lagrangian乘数$\lambda$也包括罚函数$\varepsilon_N>0$。但与前述Lagrangian乘数法和罚函数法不同的是，在该算法中
* Lagrangian乘数$\lambda$不是自变量
* 罚函数$\varepsilon_N$可以不是很大的数值
-----------------------------------------------------------------
注： 递增Lagrangian乘数法可以看成是罚函数法或Lagrangian乘数法的变种。如从Lagrangian乘数法出发，我们将（2.1）改写为
$J=1/2K x^2 -F x +\lambda (x+g_0)+(\lambda^*-\lambda) (x+g_0)$
这里，我们记$\lambda^*$为Lagrangian乘数的准确值。考虑 $\lambda^*-\lambda$为Lagrangian乘数的增分并将此式与（3.1）比较，可以看出
$\varepsilon_N (x+g_0)\approx \lambda^*-\lambda$
如此Lagrangian乘数的准确值可以如下逼近
$ \lambda^*\approx\lambda+\varepsilon_N (x+g_0)$
对应于下面的（3.5）式。
-----------------------------------------------------------------
泛函J的极值条件为其变分为零，即
(3.2) $\delta J=K x\delta x+ (\lambda+\varepsilon_Nx)\delta x -(F-\varepsilon_Ng_0) \delta x=0$
在此由于Lagrangian乘数$\lambda$不是自变量，我们不需要对其变分。基于变分$\delta x$的任意性，要求
(3.3) $(K+\varepsilon_N)x+\lambda-(F+\varepsilon_Ng_0)=0$
由此可得
(3.4) $ x=(F-\varepsilon_Ng_0-\lambda)/(K+\varepsilon_N) $
但是这里$\lambda$也是待定值，我们需要用迭代方法解此方程。
* 假设Lagrangian乘数的初值$\lambda_0=0$，n=0
* 带入上式求得$x_{n+1}$
* 判断接触条件$g_0+x>=0$是否满足。如不满足，则更新$\lambda$，
(3.5) $\lambda_{n+1}\approx \lambda_n+\varepsilon_N(g_0+x)$
令n->n+1; 并返回上一步重求x。
在上述计算过程中，Lagrangian乘数$\lambda$从零出发逐渐递增，这是我们称之为递增Lagrangian乘数法的原因。同时，x的值从$ x=(F-\varepsilon_Ng_0)/(K+\varepsilon_N) $出发逐渐递减，当$g_0+x$接近于0时计算结束。
可以看出，这里给出的x的初值与上述罚函数法给出的公式相同。因此，如果$\varepsilon_N$的值足够大，只需一步计算就可以使$g_0+x$趋于0，如此递增Lagrangian乘数法将退化为罚函数法。但是在递增Lagrangian乘数法中罚函数$\varepsilon_N$一般不是一个很大的数值，其原因将在下节讨论。
(4.1) $J=1/2K x^2 -F x -r log(g_0+x)=0$
这里$G(g)=log(x+g_0)$称为围墙（barrier）函数，因为该函数在未接触区域为一正值，但当其靠近接触面时急剧增大，趋近接触面时，其值趋向于正无穷大。它像围墙一样阻止迭代点越出约束边界。围墙函数可以有不同形式，如$G(g)=1/(g_0+x)$. 另外系数r>0的作用类似于前述罚函数。但在此它是个待定变数。在下面（4.3）式中可以看到，当r趋于零时，泛函J趋于极值。
泛函J的极值条件为其变分为零，即
(4.2) $\delta J=K x\delta x -(F+r/(g_0+x)) \delta x$
基于变分$\delta x$的任意性，要求
(4.3) $\begin{eqnarray}
Kx-\lambda-F=0 \\
\lambda（g_0+x）=r
\end{eqnarray}$
在这里我们导入了一个新变量$\lambda=r/(g+x)$. 可以看出该式与方程（2.3）非常类似，当r=0时两者等价。该方程组有三个变量x，r和$\lambda$, 需要用迭代方法求解。
* 取初值r〉0，递减系数$0<\gamma<1，\lambda_0=0, x_0=0$
* 带入(4.3)式求得x和$\lambda$. 注意到这是一个非线性方程组，需要将其线性化
(4.4) $\begin{eqnarray}
K\Delta x-\Delta \lambda=F-Kx+\lambda \\
\lambda \Delta x+(g_0+x)\Delta \lambda =r-\lambda(g_0+x)
\end{eqnarray}$
采用Newton迭代法求解。从此联立方程中可以直接求得x和$\lambda$的增量（primal-dual法），也可仅以x为自变量，靠迭代求解$\lambda$的增量（primal算法）。
* 令$r=r\gamma$, 回到上一步更新x和$\lambda$，直到r趋于0。
(3.1) $J=1/2K x^2 -F x +\lambda (x+g_0)+1/2\varepsilon_N (x+g_0)^2$
注意到该泛函既包括Lagrangian乘数$\lambda$也包括罚函数$\varepsilon_N>0$。但与前述Lagrangian乘数法和罚函数法不同的是，在该算法中
* Lagrangian乘数$\lambda$不是自变量
* 罚函数$\varepsilon_N$可以不是很大的数值
-----------------------------------------------------------------
注： 递增Lagrangian乘数法可以看成是罚函数法或Lagrangian乘数法的变种。如从Lagrangian乘数法出发，我们将（2.1）改写为
$J=1/2K x^2 -F x +\lambda (x+g_0)+(\lambda^*-\lambda) (x+g_0)$
这里，我们记$\lambda^*$为Lagrangian乘数的准确值。考虑 $\lambda^*-\lambda$为Lagrangian乘数的增分并将此式与（3.1）比较，可以看出
$\varepsilon_N (x+g_0)\approx \lambda^*-\lambda$
如此Lagrangian乘数的准确值可以如下逼近
$ \lambda^*\approx\lambda+\varepsilon_N (x+g_0)$
对应于下面的（3.5）式。
-----------------------------------------------------------------
泛函J的极值条件为其变分为零，即
(3.2) $\delta J=K x\delta x+ (\lambda+\varepsilon_Nx)\delta x -(F-\varepsilon_Ng_0) \delta x=0$
在此由于Lagrangian乘数$\lambda$不是自变量，我们不需要对其变分。基于变分$\delta x$的任意性，要求
(3.3) $(K+\varepsilon_N)x+\lambda-(F+\varepsilon_Ng_0)=0$
由此可得
(3.4) $ x=(F-\varepsilon_Ng_0-\lambda)/(K+\varepsilon_N) $
但是这里$\lambda$也是待定值，我们需要用迭代方法解此方程。
* 假设Lagrangian乘数的初值$\lambda_0=0$，n=0
* 带入上式求得$x_{n+1}$
* 判断接触条件$g_0+x>=0$是否满足。如不满足，则更新$\lambda$，
(3.5) $\lambda_{n+1}\approx \lambda_n+\varepsilon_N(g_0+x)$
令n->n+1; 并返回上一步重求x。
在上述计算过程中，Lagrangian乘数$\lambda$从零出发逐渐递增，这是我们称之为递增Lagrangian乘数法的原因。同时，x的值从$ x=(F-\varepsilon_Ng_0)/(K+\varepsilon_N) $出发逐渐递减，当$g_0+x$接近于0时计算结束。
可以看出，这里给出的x的初值与上述罚函数法给出的公式相同。因此，如果$\varepsilon_N$的值足够大，只需一步计算就可以使$g_0+x$趋于0，如此递增Lagrangian乘数法将退化为罚函数法。但是在递增Lagrangian乘数法中罚函数$\varepsilon_N$一般不是一个很大的数值，其原因将在下节讨论。
1.4 内点法
内点法是从可行域（即未接触区域g>0）内某一初始内点出发，在可行域内求极值的方法。其泛函J为(4.1) $J=1/2K x^2 -F x -r log(g_0+x)=0$
这里$G(g)=log(x+g_0)$称为围墙（barrier）函数，因为该函数在未接触区域为一正值，但当其靠近接触面时急剧增大，趋近接触面时，其值趋向于正无穷大。它像围墙一样阻止迭代点越出约束边界。围墙函数可以有不同形式，如$G(g)=1/(g_0+x)$. 另外系数r>0的作用类似于前述罚函数。但在此它是个待定变数。在下面（4.3）式中可以看到，当r趋于零时，泛函J趋于极值。
泛函J的极值条件为其变分为零，即
(4.2) $\delta J=K x\delta x -(F+r/(g_0+x)) \delta x$
基于变分$\delta x$的任意性，要求
(4.3) $\begin{eqnarray}
Kx-\lambda-F=0 \\
\lambda（g_0+x）=r
\end{eqnarray}$
在这里我们导入了一个新变量$\lambda=r/(g+x)$. 可以看出该式与方程（2.3）非常类似，当r=0时两者等价。该方程组有三个变量x，r和$\lambda$, 需要用迭代方法求解。
* 取初值r〉0，递减系数$0<\gamma<1，\lambda_0=0, x_0=0$
* 带入(4.3)式求得x和$\lambda$. 注意到这是一个非线性方程组，需要将其线性化
(4.4) $\begin{eqnarray}
K\Delta x-\Delta \lambda=F-Kx+\lambda \\
\lambda \Delta x+(g_0+x)\Delta \lambda =r-\lambda(g_0+x)
\end{eqnarray}$
采用Newton迭代法求解。从此联立方程中可以直接求得x和$\lambda$的增量（primal-dual法），也可仅以x为自变量，靠迭代求解$\lambda$的增量（primal算法）。
* 令$r=r\gamma$, 回到上一步更新x和$\lambda$，直到r趋于0。
2 算法评价
2.1 算法评价基准
评价一个数值算法的基本判据有
* 精度
* 收敛性
* 计算量
* 内存使用量
* 并行算法的相容性
对于有限元计算来说，其大部分计算时间消耗于联立方程的求解计算。求解联立方程大致可分为两类* 收敛性
* 计算量
* 内存使用量
* 并行算法的相容性
* 直接法
* 以CG法为代表的迭代法
其比较如下
直接法 迭代法
------------------------------------------------------------------- ---------
* 精度 好 可控
* 收敛性 不存在此问题 依赖于问题类型
* 计算量 较大 较小
* 内存使用量 较大 较小
* 并行算法的相容性 不好 好
从上表可以看出，对于大型计算而言，迭代法是较好的选择。在本文中将考虑各种算法对迭代法的收敛性的影响。如上所述，迭代法的收敛性依赖于问题类型，其具体的指标则是联立方程A的条件数。简单地，条件数可以理解为矩阵对角项的最大值与最小值的比
$cond(A) = max|A_{diag}|/min|A_{diag}|$
条件数增大则迭代法收敛性变差，这意味着迭代次数和计算时间的增加甚至计算发散。
另外，我们希望方程是正定的，主对角集中的，这些都可以提高迭代法求解器的收敛性。
2.2 算法评价
罚函数法:缺点： 近似解，要求极大罚函数值， 使方程组条件数变差。
优点： 计算量不大，实现最简单。
Lagrangian乘数法:
缺点： 增加了计算变量（计算量增加），方程性能变差（参见式（2.3）,可以看出其导入了零对角项，该方程变为非正定方程）。
优点： 精确解。
递增Lagrangian乘数法:
缺点： 需要迭代求解Lagrangian乘数（计算量增加）。
优点： 回避了计算变量的增加；由于选用较小的罚函数值，可以回避或减缓方程组条件数的恶化；对于接触问题解析来说，可以利用此方法把非对称的接触刚性矩阵对称化，大幅度节省内存和计算时间。
内点法:
缺点： 不能处理初始点已经位于接触面内的问题；采用primal-dual算法增加了计算变量；需要迭代求解。
优点： 内点法应用于接触问题解析主要是其可以回避解析过程中，接触点—〉非接触-〉接触，即所谓接触active set的变化，因为该方法对接触面附近的所有点都加上了约束。因此可以期待其提高计算的收敛性。
3 递增Lagrangian乘数法的实装方法（待续）
*1 Augmented Lagrangian的译法有多种。因为 augment的原意为 To make (something already developed or well under way) greater, as in size, extent, or quantity （http://www.thefreedictionary.com/augment），Make (something) greater by adding to it; increase (http://www.oxforddictionaries.com/definition/english/augment). 而Augmented Lagrangian法中，Lagrangian乘数的值正是逐渐增大的，因此提出该译法。
*2 接触条件一般在接触点为零点，接触面法线方向为正向的局部坐标系下表述。由于本例题只有一个自由度，本文简单地取局部坐标系为全局坐标。
*3 由于习惯的不同，(3.1)式也可以写为$J=1/2K x^2 -F x -\lambda (x+g_0)+1/2\varepsilon_N (x+g_0)^2$。这时得到的$\lambda$与(3.1)式相比差一个负号。
Compile p_arpack and arpack-ng in mingw
1 arpack-ng
The configure script of arpack-ng is not that smart. You need define several settings by hand like
export MPIF77=gfortran
export FFLAGS=-I/c/mpich2/include
./configure --with-blas=/c/libs/libblas.a --with-lapack=/c/libs/liblapack.a --enable-mpi
where absolute location of blas and lapack is needed. Besides, the MPICH2 library libfmpich2g.a is not considered in the configure script, you need add something like
LIBS="-lfmpich2g $LIBS"
cat > conftest.$ac_ext <<_ACEOF
program main
call MPI_Init
end
_ACEOF
if ac_fn_f77_try_link "$LINENO"; then :
ac_cv_lib_fmpich_MPI_Init=yes
else
ac_cv_lib_fmpich_MPI_Init=no
fi
to check the existence of MPILIBS
2. Pararpack
The "mpif.h" file in the pararpack pakage is too outdated to be used. You must delete this files and modify ARmake.inc to indicate MPI settings. Like
MPIDIR = /c/MPICH2
MPILIBDIR = $(MPIDIR)/lib
MPILIBS = -L$(MPILIBDIR) -lfmpich2g -lmpi
FFLAGS = -O -I/c/MPICH2/include
The configure script of arpack-ng is not that smart. You need define several settings by hand like
export MPIF77=gfortran
export FFLAGS=-I/c/mpich2/include
./configure --with-blas=/c/libs/libblas.a --with-lapack=/c/libs/liblapack.a --enable-mpi
where absolute location of blas and lapack is needed. Besides, the MPICH2 library libfmpich2g.a is not considered in the configure script, you need add something like
LIBS="-lfmpich2g $LIBS"
cat > conftest.$ac_ext <<_ACEOF
program main
call MPI_Init
end
_ACEOF
if ac_fn_f77_try_link "$LINENO"; then :
ac_cv_lib_fmpich_MPI_Init=yes
else
ac_cv_lib_fmpich_MPI_Init=no
fi
to check the existence of MPILIBS
2. Pararpack
The "mpif.h" file in the pararpack pakage is too outdated to be used. You must delete this files and modify ARmake.inc to indicate MPI settings. Like
MPIDIR = /c/MPICH2
MPILIBDIR = $(MPIDIR)/lib
MPILIBS = -L$(MPILIBDIR) -lfmpich2g -lmpi
FFLAGS = -O -I/c/MPICH2/include
订阅：
博文 (Atom)
VS Code下cmake, c++编译，调试环境的构成步骤
1 下载必须extension 按[Ctrl+Shift+X]打开extension窗口，选择安装"C/C++", "CMake", "CMake Tools" 2 在VSCode下打开作业目录 ...
-
学过连续体力学(包括固体力学,流体力学等)都知道应力张量是对称的(下面将说明它是有条件的)。 对于初学者， 这似乎有些难于理解， 如 http://forum.simwe.com/forum.php?mod=viewthread&tid=503234&hig... -
1. 什么是体积锁死现象 图1 平面应变模型 在图1所示的又两个三角形单元组成的平面应变问题中，如果变形体材料是不可压缩的，那么我们不管在加载点施加多大的力，从有限元法计算得到的所有节点的所有位移都为零。因为任意位移都会导致 两个三角形中的一个的体积... -
假设时间步i中的所有物理量X i 已知, 如果时间步i+1中的所有物理量X i+1 可以表为X i 的显函数 （1） X i+1 = F(X i ) （注；本文以下标i，i+1等表示时间步） 则该物理量X可以显式计算。而如果 （2） ...