上善若水: 2014

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

有限元接触解析方法（之一）

罚函数法，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 $
这是该问题的精确解。

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$一般不是一个很大的数值，其原因将在下节讨论。

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

内存误用的检查方法

检查是否有内存遗漏的简单方法
i = system( "wmic os get FreePhysicalMemory /value" )
该命令显示当前可用物理内存。

Linux下可以用valgrind去特定内存遗漏
valgrind --error-limit=no --leak-check=full --show-reachable=no myprog arg1 arg2 2>&1 | tee valgrind.log
详细可参见http://valgrind.org/docs/manual/mc-manual.html

Windows下可用DrMemory : http://www.drmemory.org/

其他
Wine and Valgrind http://wiki.winehq.org/Wine_and_Valgrind
Visual Leak Detector for Visual C http://vld.codeplex.com/
Very sleepy http://www.codersnotes.com/sleepy
VMMap http://technet.microsoft.com/en-us/sysinternals/dd535533
D.U.M.A. http://duma.sourceforge.net/

梯度是向量还是1形式？

一直以为梯度是向量，可最近看Schutz的A first course in general relativity, 该书的3.3节下有一小节的标题是Gradient of a function is a one-form称梯度是1形式。

$df=\frac{\partial f}{\partial x} dx+\frac{\partial f}{\partial y} dy+\frac{\partial f}{\partial z} dz$

Schutz直接定义由成份$\partial f/\partial x，\partial f/\partial y， \partial f/\partial z$构成的量为梯度，这是1形式。

但是，一般情况下，梯度定义为

$\nabla f=g^{ik} \frac{\partial f}{\partial x^k} \frac{\partial}{ \partial x^i}$

其成分为$g^{ik}\frac{\partial f}{\partial x^k}$,这是向量.

参考文献
1. B. Schutz: A first course in general relativity, Cambridge Univ, 2009
2. http://en.wikipedia.org/wiki/Gradient

Lagrangian to Hamiltonian = TQ to T*Q = Riemann to Sympletic

假设系统有Lagrangian
$L=L(q^1,...,q^n,\dot{q}^1,...\dot{q}^n,t)$
并满足 Euler-Lagrange 方程
(1)    $d/dt \partial L/\partial\dot{q}^i-\partial L/\partial q^i=0$
可以看出，上述微分方程处于2n维空间TQ，其坐标为
      $q^1,...,q^n,\dot{q}^1,...\dot{q}^n$
可以解释为构形空间Q的切丛。
通常设势能V为$q^i$的函数，动能$T=1/2g_{ij}\dot{q}^i\dot{q}^j$。这里$g_{ij}$为Riemann测度张量。因此方程（1）为Riemann空间中的极值问题，特别是V=0时得到的就是Riemann空间中的测地线。

    现在使用Legendre变换从切丛TQ转换至余切丛 T*Q
$\omega$: TQ --> T*Q
(q,$\dot{q}$)--> (q,p)
这里 $p_i=\partial L/\partial\dot{q}^i, i=1,...,n.$. 这时由Lagrangian函数描述的处于切丛 TQ的系统将转换为由下述 Hamiltonian函数 H描述的处于相空间 T*Q 的系统
       $H(p,q,t):=p_i\dot{q}^i-L(q,\dot{q},t) \; with\; p_i=\partial L/\partial\dot{q}^i,$
Lagrange方程 (1)也转换为 Hamilton方程
(2)          $\dot{q}^i=\partial H/\partial p_i, \dot{p}_i=-\partial H/\partial q^i$
写成 Poisson 括号形式则为
(3)         $\dot{q}^i=\{q^i,H\}, \dot{p}_i=\{p_i,H\}$

积分方程（2),我们可以得到p,q在相空间的流动。进一步，我们研究相空间的Sympletic结构。
考虑Q上的可微函数H，其微分为一1形式$dH=\partial H/\partial q^idq^i+\partial H/\partial p_idp_i$.将其拉回到TQ则得到Hamiltonian矢量$V_H=\omega^{-1}(dH)$,写成内积(interior product)则为
$\omega(V_H,*)=i_{V_H}\omega=dH$
由于我们要求变换$\omega: TQ-->T^*Q$一一对应，变换$\omega$必须非退化。另外我们要求Hamiltonian守恒，即
$dH(V_H)=\omega(V_H,V_H)=0$
则$\omega$须为2-形式，即$\omega$具有Sympletic结构（非退化的封闭的2-形式）
$\omega=dp_i\wedge dq^i$

References:
1. R. Berndt: An Introduction to Sympletic Geometry, American Mathematical Society, 2001
2. Topics in Representation Theory:Hamiltonian Mechanics and Symplectic Geometry
http://www.math.columbia.edu/~woit/notes22.pdf
3. SYMPLECTIC VECTOR FIELDS http://www-personal.umich.edu/~wangzuoq/437W13/Notes/Lec%2034.pdf

QWidget: Must construct a QApplication before a QWidget

There are instances that you download official examples from http://www.vtk.org/Wiki/VTK/Examples/Cxx/Qt. but got error like "QWidget: Must construct a QApplication before a QWidget" during runtime. It is most probably due to the debug and release conflict between Qt libaray and your exe file. Just rebuild your exe file by set CMAKE_BUILD_TYPE to Debug or Release correspondingly. It can be done by modifying CMakefile.txt as
--------------------------------------------
set(CMAKE_BUILD_TYPE Debug) or set(CMAKE_BUILD_TYPE Release)
-------------------------------------------
or running cmake like
-------------------------------------------
cmake -DCMAKE_BUILD_TYPE=Release
-------------------------------------------

应力张量为什么是对称的? 为什么是不对称的?(其二)

4 什么应力张量是不对称的
   如前文所指出， Cauchy应力张量的对称性得至于角动量平衡方程(2). 如果方程(2)不成立，如将其改写为
$\int_\Omega x\times\rho \dot{u}dv=\int_\Omega x\times b dv+\int_{\partial\Omega} x\times t ds + \int_\Omega m dv$                     (10)
则Cauchy应力张量就会变得不对称性。这里增加的一项m的物理意义为外部体力矩。这是一种假设，至于如何在尺寸为零的点上施加力矩这样的问题比较难于回答，在这里我们先放过了。
   上述在古典连续体力学导入新的变量的想法其实并不新鲜，它最早出现在1909年（[17]）。Cosserat理论可以说是现在渐渐变得时髦的nonlocal continuum或 generalize continuum mechanics的鼻祖。在上世纪60年代对类似理论有详细的研究。这样构筑的力学体系于古典连续体力学大不相同，它们有自己的平衡方程，本构方程，Lagrangian或Hamilton. 在这些理论体系中，Cauchy应力有些是对称的（如couple stress theory, strain graident theory），有些是不对称得（如micropolar theory)。本文无意讨论这些问题，下述文献[18]-[20]叙述的理论中的Cauchy应力都是不对称的。

5 结言：模型与现实
   在本文中，我们提到了至少两类连续体模型：一类的Cauchy应力是对称的，另一类则是不对称的。我们在构筑模型时，第一要求模型体系是自洽的，即不是自相矛盾的。经典连续体力学就是一个这样的严格的数学体系，在该体系下不能也不可能得到不对称的Cauchy应力。第二我们要求模型可以解释已知并且可以预测物理现象。进一步，如果模型的结果与现实不相符，该模型就需要修正或者放弃重构。由于经典连续体力学在某些方面的预测结果的精度不够，nonlocal continuum或 generalize continuum mechanics正是对经典连续体力学的修正。

6 题外的话：为什么应变张量是对称的？
在（8）式中，我们看到应变能为Cauchy应力和应变速率$\varepsilon=gradu$的积。这里定义的应变速率不是对称的,它可以分解为对称部分和反对称部分的和。考虑到Cauchy应力的对称性, 应力张量与反对称张量的乘积为零，应变能实际上只和$gradu$的对称部分有关，即$\sigma:gradu=\sigma:sym(gradu)$。这样应变速率可以定义为$\varepsilon=sym(gradu)$. 这是常见的应变张量定义方法。

参考文献

17 Cosserat; E. & F., Théorie de Corps déformables, Paris, 1909
18 Eringen, A.C; Theory of micropolar elasticity, Defense Technical Information Center, 1967
19 Nowacki, W; The linear theory of micropolar elasticity, 1972
http://bcpw.bg.pw.edu.pl/Content/970/Micropolar_Elasticity_str1_43.pdf
20 Nowacki, W; Theory of Asymmetric Elasticity, Elsevier Science & Technology, 1986

Build scotch6.0 in mingw64

1. Prerequisite
    Following tools and libraries install in mingw64: flex, yacc, zlib, pthread, mpi

2. Modify the Makefile.inc
    Based upon file Makfile.inc.i686_pc_mingw32, modify in follows
1)  Add -Drestrict=__restrict into CFLAGS_DEF
2) Modify -lmpich2mt in LDFLAGS into -lmpicxx
3) Modify pthread library into what you have installed

3. Add following into common.h
#ifndef HAVE_STDINT_H
#define HAVE_STDINT_H   1
#endif /* HAVE_STDINT_H */
#ifndef HAVE_UINT_T
#define HAVE_UINT_T   1
#endif /* HAVE_UINT_T */

Then, everything ok!

应力张量为什么是对称的? 为什么是不对称的?（其一）

学过连续体力学(包括固体力学,流体力学等)都知道应力张量是对称的(下面将说明它是有条件的)。对于初学者，这似乎有些难于理解，如
http://forum.simwe.com/forum.php?mod=viewthread&tid=503234&highlight=%E5%89%AA%E5%BA%94%E5%8A%9B%E4%BA%92%E7%AD%89

奇妙的是，有些专家也对此持疑。如下图的例子([1] p.110)就出自中国力学院士的专著，他给出了几乎与上面连接给出的同样的图来说明应力是不对称的。这就有点哭笑不得了。

本文试图对应力张量的对称性做一个较直观的但有些深度的解释。在此之前，先指出上图论证的问题- 我们不能对连续体的一点（它的体积→零，表面积→为零）的某一面施加外力，也就是说不能直接指定应力分量。外力边界条件，一般有称为第二边界条件（[2]，p50），在这里你可以给出外力矢量（三个分量），而不是六个应力分量。因此上述论理是不成立的。

1 什么是应力？
   要理解应力当然先要理解什么是力。
  很多人认为牛顿第二定律给出了力的定义。但牛顿第二定律给出的力的计算式has no independent meaning（[4]）。力的定义有一定的任意性，它也许毫无道理，但并不奇怪（It may be gratuitous, but it is not bizarre[4]）。Feynman先生认为上述定义是无用的（[3]&12.1:The Newtonian statement above, however, seems to be a most precise definition of force, and one that appeals to the mathematician; nevertheless, it is completely useless），似乎不存在力的精确定义（[3]&12.1: If you insist upon a precise definition of force, you will never get it!)。既然两位诺贝尔奖获得者都这么说。我们还是放弃在这里定义力，假设它是一种基本物理量为好。但是要注明的是，在现代物理学中，力不是基本物理量，它一般被理解为能量的空间导数或动量的时间导数。此时应力的定义也相应解释为能量密度的导数等(如[5])。但本文不采用这种不直观的定义方法。
   记过连续体一点x的任意切面（法线方向n）的表面力矢量场为f（n，x）。Cauchy定理指出
$f（n，x）=\sigma(x)n$
这里的$\sigma(x)$即为Cauchy应力张量。

几乎所有的连续体力学教程都会写下上述Cauchy定理的证明，但严格的少见。建议阅读文献[6]p101-105，或[7]p26-27, [8]-[9]。Cauchy定理仍然有议论的余地，如放松定理成立的光滑条件，考虑上述表面力矢量场是法线方向的空间导数的函数，将其导入微分几何学等（[10]-[14]）。

2 为什么Cauchy应力张量是对称的
考虑一物体，其动量为$\int_\Omega \rho udv$, 所受体积力为$\int_\Omega b dv$ ，面力为$\int_{\partial\Omega} t ds$。则由牛顿第二定律

$\int_\Omega \rho \dot{u}dv=\int_\Omega b dv+\int_{\partial\Omega} t ds$ (1)

其角动量平衡方程为

$\int_\Omega x\times\rho \dot{u}dv=\int_\Omega x\times b dv+\int_{\partial\Omega} x\times t ds$ (2)

由散度定理$\int_{\partial\Omega} t ds=\int_\Omega div\sigma dv$, 将其代入方程（1）得

$\int_\Omega (div\sigma+b-\rho \dot{v})dv=0$ (3)

由于该式在连续体内任意一点都必须成立，得到平衡方程

$div\sigma+b-\rho \dot{v}=0$ (4)

   从方程（2）则可以得到应力张量的对称性。方程（2）的最后一项$\int_{\partial\Omega} x\times t ds=\int_{\partial\Omega} x\times\sigma n ds=\int_\Omega(x \times div\sigma+\epsilon:\sigma^T)dv$     (5)
将此式代入（2）并使用方程（4），即可得到
$\epsilon:\sigma^T=0$             （6）
即Cauchy应力张量是对称。
如上所述，Cauchy应力张量的对称性来源于角动量平衡条件，如果Cauchy定理成立，则Cauchy应力必然是对称的。

3 单独定义的应力没有什么实用意义（has no independent meaning）
   上节的平衡方程（4）和Cauchy应力张量的对称性条件实际上只是牛顿定律应用于连续体时的再述。上述方程并未给出关于连续体变形的任何信息。为此，我们需要导出于Cauchy应力张量
共轭的量，该张量与Cauchy应力张量的积表征连续体的变形能。或用现代物理学的语言来说，可以用来构筑变形体系统的Lagrangian或Hamilton。
与上节相似，体积力所作的功为$\int_\Omega b．vdv$。表面力所作的功为$\int_{\partial\Omega} t．uds$。则外力的总功W为
$W=\int_\Omega b．udv+\int_{\partial\Omega} t．udv$       (7)
其中，表面力的功$\int_{\partial\Omega} t．udv=\int_{\partial\Omega} (\sigma n)．udv=\int_\Omega div(\sigma^Tu)dv=\int_\Omega (div\sigma\cdot u+\sigma:gradu)dv$。因此
$W=\int_\Omega (\rho\dot{u}\cdot{u}+\sigma:gradu)dv+\int_{\Omega}(div\sigma+b-\rho \dot{u})\cdot uds=\int_\Omega (\rho\dot{u}\cdot{u}+\sigma:gradu)dv$    (8)
该式的倒数第二项为变形体的动能，倒数第一项为变形能。记$\varepsilon=gradu$为应变速率，该应变速率张量为Cauchy应力的共轭量。

3.1 应力应变张量在不同构型(configuration)下的表达
    由于连续体的形状是变化的，如同物理量可以在不同的坐标系下表示一样，应力，应变也可以在不同构型下表达。下面是一个例子
    考虑现在构型$\Omega$ 下的变形能，它由Cauchy应力，应变速率表达. 下面我们将其变换到构型$\Omega_0$. 两构型间的两点变换张量为F, Jocabian为J。
$W=\int_\Omega\sigma:\varepsilon dv=\int_\Omega\sigma:（\dot F F^{-1})dv=\int_\Omega\sigma F^{-T}:\dot Fdv=\int_{\Omega_0}J\sigma F^{-T}:\dot FdV$
   一般定义$P=J\sigma F^{-T}$为第一种Piola-Kirchhoff应力。它与$\dot F$共轭。另外由于F不是对称的，所以第一种Piola-Kirchhoff应力不对称。

3.2 从能量表达式（7）推出平衡方程
    从现代物理学的观点来看，整个连续体力学都可以建立于系统的能量表达式。我们可以从能量表达式（7）出发，利用标架不变(frame-indifference)得出各守恒定理如式（2），（3）。可参见如[15]，[16]）。也许对基础坚实的古典连续体力学来说这样做的意义不大，但如要将古典连续体力学加以推广，这是一个可靠的工具。

参考文献
1. 陈至达: 理性力学, 2000，重庆出版社
2 穆什海里什维利: 数学弹性力学的几个基本问题
3 费曼物理学讲义， http://www.feynmanlectures.caltech.edu/
4 Frank Wilczek: Whence the force of F=ma? http://ctpweb.lns.mit.edu/physics_today/phystoday/%20Whence_cshock.pdf
5 Robert G. Brown, 2013: Symmetric stress tensor; http://www.phy.duke.edu/~rgb/Class/Electrodynamics/Electrodynamics/node147.html
6. M.E.Gurtin: A introduction to Continuum mechanics, 1981, Academic press
7 J.T.Oden:A short course on nonlinear continuum mechanics, 2008 http://users.ices.utexas.edu/~arbogast/cam397/oden0908.pdf
8 Miroslav Šilhavy：On Cauchy's stress theory, http://www.bdim.eu/item?fmt=pdf&id=RLIN_1990_9_1_3_259_0
9 R.L.Fosdick, E.G.Virga: A viariaional proof of stress theroem of Cauchy, Archives of Rational mechanica and analysis, 1998, p95-103
10 G. RODNAY AND R. SEGEV: Cauchy's flux theorem in light of geometric integration theory. http://www.bgu.ac.il/~rsegev/Papers/FluxGeomIntegration.pdf
11 Francesco dell’Isola et al; How contact interactions may depend on the shape of Cauchy cuts in N-th gradient continua: approach; http://hal.archives-ouvertes.fr/docs/00/66/23/76/PDF/dellisola_seppecher_madeo.pdf
12 W. Noll: Thoughts on the concept of stress. http://repository.cmu.edu/cgi/viewcontent.cgi?article=1015&context=math
13 E.Kanso等:On geometric character of stress in continuum mechanics http://upcommons.upc.edu/e-prints/bitstream/2117/8516/1/kanso_on-the-geometric_2007.pdf
14 C.A.Trusdell: Cauchy and the modern mechanics of continua
http://www.persee.fr/web/revues/home/prescript/article/rhs_0151-4105_1992_num_45_1_4229
15 P.Germain; The method of virtual power in continuum mechanics, Part II; SIAM J. Appl. Math, Vol. 25(1973), p556-575
16 G D Piero; On the method of virtual power in coninuum mechanics. J. Mech. Mater. Struct., VOl.4(2009), p281-292 http://msp.org/jomms/2009/4-2/jomms-v4-n2-p07-p.pdf
17   G D Piero; Virtual power, pseudovalance and the law of action and reaction； http://www.fyffm2010.cnrs-mrs.fr/PDFs/Del_Piero_Gianpietro.pdf

订阅：博文 (Atom)

上善若水

Build Opensees in Visual studio 2012

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

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