Fortran Coder

标题: 坐标变换程序运行得不到结果 [打印本页]

作者: hipeilei 时间: 2015-11-12 10:10
标题: 坐标变换程序运行得不到结果
Fortran菜鸟一枚，写了一个坐标变换的Fortran程序，大致的过程就是随便给出一个整体坐标系下的一个点（x，y），通过计算得到该点在新的局部坐标系下的对应点（m，n），局部坐标系就是四边形等参变换的局部坐标系，各点坐标如下图右边所示。算法采用的等参变换的原理，如下图，整个过程用了雅可比矩阵，使用迭代的方法去判断，直到满足精度要求就认为找到相应局部坐标系下对应的点，这里选择的整体坐标下的四边形各点坐标为x1=-1,y1=-1,x2=1,y2=-1,x3=1,y3=1,x4=-1,y4=1，要判断的点假设为x=0,y=0。这样一来，如果程序没错就会给出m=0，n=0的结果。程序编译运行都没问题，但运行以后发现不论x，y取何值，程序输出结果始终都是和初值m=0，n=0相等。现贴出代码，请各位大神指点一二！

[Fortran] 纯文本查看 复制代码

program main

implicit none

real::x1=-1,y1=-1,x2=1,y2=-1,x3=1,y3=1,x4=-1,y4=1

real:: x=0,y=0

real::m=0,n=0

!read(*,*)x1,y1,x2,y2,x3,y3,x4,y4

!read(*,*)x,y

call iteration(x1,y1,x2,y2,x3,y3,x4,y4,x,y,m,n)

write(*,*)m,n

pause

end



subroutine Jacobi(x1,y1,x2,y2,x3,y3,x4,y4,m,n,a,b,c,d)

implicit none

real:: a,b,c,d,m,n

real:: x1,y1,x2,y2,x3,y3,x4,y4

a=(x1*(-1+n)+x2*(1-n)+x3*(1+n)+x4*(-1-n))/4

b=(y1*(-1+n)+y2*(1-n)+y3*(1+n)+y4*(-1-n))/4

c=(x1*(-1+m)+x2*(-1-m)+x3*(1+m)+x4*(1-m))/4

d=(y1*(-1+m)+y2*(-1-m)+y3*(1+m)+y4*(1-m))/4

end subroutine Jacobi



subroutine shapefunction(m,n,N1,N2,N3,N4)

implicit none

real:: m,n,N1,N2,N3,N4

N1=(1-m)*(1-n)/4

N2=(1+m)*(1-n)/4

N3=(1+m)*(1+n)/4

N4=(1-m)*(1+n)/4

end subroutine shapefunction



subroutine iteration(x1,y1,x2,y2,x3,y3,x4,y4,x,y,m,n)

real::x0,y0,condition1,condition2,m0,n0,m,n

real::N1=0,N2=0,N3=0,N4=0

real::a=0,b=0,c=0,d=0

real,dimension(2,2)::J

do 

  call shapefunction(m,n,N1,N2,N3,N4)

  x0=N1*x1+N2*x2+N3*x3+N4*x4

  y0=N1*y1+N2*y2+N3*y3+N4*y4

  call Jacobi(x1,y1,x2,y2,x3,y3,x4,y4,m,n,a,b,c,d)

  J(1,1)=d/(a*d-b*c)

  J(1,2)=-b/(a*d-b*c)

  J(2,1)=-c/(a*d-b*c)

  J(2,2)=a/(a*d-c*b)

  m0=m+J(1,1)*(x-x0)+J(1,2)*(y-y0)

  n0=n+J(2,1)*(x-x0)+J(2,2)*(y-y0)

  condition1=m0-m

  condition2=n0-n

   if(abs(condition1)<0.001.and.abs(condition2)<0.001) then

     m=m0

     n=n0

     exit

  else

     m=m0

     n=n0

  endif

enddo

end subroutine iteration

作者: fcode 时间: 2015-11-12 10:54
既然是迭代，那m n 应该有初值。

作者: hipeilei 时间: 2015-11-12 11:20
本帖最后由 hipeilei 于 2015-11-12 11:25 编辑

fcode 发表于 2015-11-12 10:54
既然是迭代，那m n 应该有初值。

老大，这里m=o，n=0，结果是0，0，不管怎么取值永远都是0，0.

作者: li913 时间: 2015-11-12 11:35
程序运行结果不对，这种情况是逻辑错误，也就是算法写错了，只有熟悉算法才能解决。所以，详细描述你的算法才行。我百度了一下“等参变换”，还是不明白。

作者: fcode 时间: 2015-11-12 11:43

hipeilei 发表于 2015-11-12 11:20
老大，这里m=o，n=0，结果是0，0，不管怎么取值永远都是0，0.

你作弊。之前你给的代码里并没有对 m n 赋予初值。

你后来修改了。我使用你修改后的代码，得到了 0.0 和 0.0 的结果。

QQ截图20151112114541.png (25.37 KB, 下载次数: 370)

作者: hipeilei 时间: 2015-11-12 15:17

fcode 发表于 2015-11-12 11:43
你作弊。之前你给的代码里并没有对 m n 赋予初值。

你后来修改了。我使用你修改后的代码，得到了 0.0 和 ...

老大，不好意思

我当时调的时候改来改去删掉m，n初值以后忘了加上去了，我加上去以后发现是可以输出，但是不是我想要的结果啊，因为不管我的x，y如何取值，最终的输出结果都只是0，0 也就是说我的初值m，n取的是几输出的就是几，完全没有运行我中间的子程序啊

作者: hipeilei 时间: 2015-11-12 16:29

li913 发表于 2015-11-12 11:35
程序运行结果不对，这种情况是逻辑错误，也就是算法写错了，只有熟悉算法才能解决。所以，详细描述你的算法 ...

算法说起来有点长，且听我一一道来，论坛里面格式不好编辑，还是看图吧，看下图：

作者: kerb 时间: 2015-11-12 18:48
好像是有限元等参变换，似乎你搞反了，实际上是[latex]P_i(x_i,y_y)(i=1,2,3,4)[/latex]变换到[latex]\Omega_i(\xi_i,\eta_i)[/latex]空间，应该是[latex](\xi_i,\eta_i)=(\pm1,\pm1)[/latex]四种组合，但是你写成x1=,y1=...x4=,y4=...，按照你的程序，应该是m,n=-1,+1...你看看下图有助于理解吗

360截图20151112183841171.jpg (24.32 KB, 下载次数: 423)

作者: hipeilei 时间: 2015-11-12 23:15
本帖最后由 hipeilei 于 2015-11-12 23:16 编辑

kerb 发表于 2015-11-12 18:48
好像是有限元等参变换，似乎你搞反了，实际上是变换到空间，应该是四种组合，但是你写成x1=,y1=...x4=,y4=. ...

感谢回答！这就是有限元的等参变换问题，你可能误解我的意思了，也怪我没说清楚，我这里选的坐标点（如下图）是整体坐标的，就是那个不规则四边形的四个点，我选择和等参单元一样的坐标是为了便于验证结果；

有限元的的等参变换是把不规则的四边形四个节点变成规则的四个节点坐标（如下图）

规则的四边形的四个节点如图，就和你说的一样。
问题来了：（这也是我程序实现要的目的）：我这里想要做的事就是，给出一个不规则四边形里面任意一个点（上图左），找到在规则的四边形里面对应的点。过程用到了迭代的方法，这里不规则四边形的四个坐标已知，四边形内的坐标点(x0，y0)也已知，任取一个规则四边形内的坐标点（m，n）【相当于下图的kesi，yita】，我这里取的是（0，0),通过下面的式子

根据迭代的方法，知道满足一下条件：

就认为找到了对应的（kesi，yita），可是这一算法一直算不出来正确的结果，我已经修改了很多次了，这是最近一次修改，麻烦给看看：
program main
implicit none
real::x1=-1,y1=-1,x2=1,y2=-1,x3=1,y3=1,x4=1,y4=1
real:: x=0.3,y=0.3
real::m=0.2,n=0.3
real,dimension(2,2)::J
real::x0,y0,condition1,condition2,m0,n0
real:: N1=0,N2=0,N3=0,N4=0,a,b,c,d
!read(*,*)x1,y1,x2,y2,x3,y3,x4,y4
!read(*,*)x,y

call shapefunction(m,n,N1,N2,N3,N4)
  x0=N1*x1+N2*x2+N3*x3+N4*x4
  y0=N1*y1+N2*y2+N3*y3+N4*y4
  !write(*,*)x0,y0
call Jacobi(x1,y1,x2,y2,x3,y3,x4,y4,m,n,a,b,c,d)
  m0=m+J(1,1)*(x-x0)+J(1,2)*(y-y0)
  n0=n+J(2,1)*(x-x0)+J(2,2)*(y-y0)
  !write(*,*)m0,n0
  condition1=m0-m
  condition2=n0-n
  do while((condition1)<0.001.and.abs(condition2)<0.001)
   m=m0
   n=n0
  enddo
!write(*,*)a,b,c,d
!write(*,*)J(1,1),J(1,2),J(2,1),J(2,2)
write(*,*)m,n
  pause
end

subroutine shapefunction(m,n,N1,N2,N3,N4)
implicit none
real:: m,n,N1,N2,N3,N4
N1=(1-m)*(1-n)/4
N2=(1+m)*(1-n)/4
N3=(1+m)*(1+n)/4
N4=(1-m)*(1+n)/4
write(*,*),"***********This is shapefuction!*******************"
write(*,"(4f10.2)")N1,N2,N3,N4
end subroutine shapefunction

subroutine Jacobi(x1,y1,x2,y2,x3,y3,x4,y4,m,n,a,b,c,d)
implicit none
real:: a,b,c,d,m,n
real,dimension(2,2)::J
real:: x1,y1,x2,y2,x3,y3,x4,y4
a=(x1*(-1+n)+x2*(1-n)+x3*(1+n)+x4*(-1-n))/4
b=(y1*(-1+n)+y2*(1-n)+y3*(1+n)+y4*(-1-n))/4
c=(x1*(-1+m)+x2*(-1-m)+x3*(1+m)+x4*(1-m))/4
d=(y1*(-1+m)+y2*(-1-m)+y3*(1+m)+y4*(1-m))/4
  J(1,1)=d/(a*d-b*c)
  J(1,2)=-b/(a*d-b*c)
  J(2,1)=-c/(a*d-b*c)
  J(2,2)=a/(a*d-c*b)
write(*,*),"***********This is Jacobi matrix!*******************"
write(*,"(4f10.2)")J(1,1),J(1,2),J(2,1),J(2,2)
end subroutine Jacobi

作者: kerb 时间: 2015-11-12 23:42
本帖最后由 kerb 于 2015-11-12 23:46 编辑

好像可以解出两个坐标系变换，稍微麻烦一点
[latex]x=a_x+b_x\xi+c_x\eta+d_x\xi\eta[/latex]

[latex]y=a_y+b_y\xi+c_y\eta+d_y\xi\eta[/latex]
其中根据对应关系可以求出[latex]a_x ,a_y \cdots\cdots d_x, d_y[/latex]

先消去[latex]\xi\eta[/latex]这一项，然后代换，就可以化成一个关于[latex]\xi[/latex]或者[latex]\eta[/latex]的一元二次方程

作者: kerb 时间: 2015-11-13 00:16
本帖最后由 kerb 于 2015-11-13 00:22 编辑

另外用迭代法求解非线性方程组应该是这样
[latex]F(\xi,\eta)=\left\{\begin{array}{c}f_1(\xi,\eta)=a_x+b_x\xi+c_x\eta+d_x\xi\eta-x\\f_2(\xi,\eta)=a_y+b_y\xi+c_y\eta+d_y\xi\eta-y\end{array}\right.[/latex]

[latex]F(0+\delta)=F(0)+\nabla F\cdot\delta+O(\delta^2),\delta=(\Delta\xi,\Delta\eta)[/latex]所以那个雅克比矩阵是对[latex]\xi,\eta[/latex]求导

也就是：
[latex]J=\left(\begin{array}{cc}b_x+d_x\eta & c_x+d_x\xi\\b_y+d_y\eta & c_y+d_y\xi\end{array}\right)[/latex]

你检查一下是否搞混了

作者: hipeilei 时间: 2015-11-13 23:10
本帖最后由 hipeilei 于 2015-11-14 00:12 编辑

kerb 发表于 2015-11-13 00:16
另外用迭代法求解非线性方程组应该是这样

多谢指点，今天考了一天的试，到现在才看到，感谢关注。我感觉也是雅可比矩阵这里出了问题，可就是不知道到底错在什么地方了，我用maple检查了一下我的雅可比矩阵，和我代码里面的是一样的，感觉没求错。下面是maple求导的结果：

这是代码里面的式子，其中m相当于kesi，n相当于yita，a表示x对kesi求偏导，b表示y对kesi求偏导，c表示x对yita求偏导，d表示y对yita求偏导：
subroutine Jacobi(x1,y1,x2,y2,x3,y3,x4,y4,m,n,a,b,c,d)
implicit none
real:: a,b,c,d,m,n
real,dimension(2,2)::J
real:: x1,y1,x2,y2,x3,y3,x4,y4
a=(x1*(-1+n)+x2*(1-n)+x3*(1+n)+x4*(-1-n))/4
b=(y1*(-1+n)+y2*(1-n)+y3*(1+n)+y4*(-1-n))/4
c=(x1*(-1+m)+x2*(-1-m)+x3*(1+m)+x4*(1-m))/4
d=(y1*(-1+m)+y2*(-1-m)+y3*(1+m)+y4*(1-m))/4
  J(1,1)=d/(a*d-b*c)
  J(1,2)=-b/(a*d-b*c)
  J(2,1)=-c/(a*d-b*c)
  J(2,2)=a/(a*d-c*b)
write(*,*),"***********This is Jacobi matrix!*******************"
write(*,"(4f10.2)")J(1,1),J(1,2),J(2,1),J(2,2)
end subroutine Jacobi
调试的时候发现了一个新问题，就是雅可比矩阵求出来居然是0；那就说明是我的雅可比矩阵求的有问题

答主应该也是力学的吧，方便qq交流吗？我的qq：913111805

作者: hipeilei 时间: 2015-11-13 23:37

kerb 发表于 2015-11-12 23:42
好像可以解出两个坐标系变换，稍微麻烦一点

这样解出来可能会有两个根，取舍问题不好判断，有可能两个点都在规则单元里。

作者: pasuka 时间: 2015-11-15 16:50

hipeilei 发表于 2015-11-13 23:37
这样解出来可能会有两个根，取舍问题不好判断，有可能两个点都在规则单元里。 ...

碰上二义性吧？一般四边形处理云图或等值线会遇到，解决办法有很多
最直观的办法就是：四边形加条对角线切成2个三角形，假定点q不在四条边上，个么q要么在两个三角形内部，要么在对角线上，等参变换后亦如此
至于如何判断q是否在三角形内部，有个帖子已经介绍过了

作者: kerb 时间: 2015-11-18 15:50
你手头有maple，直接解析求出[latex](\xi,\eta)[/latex]的表达式，然后稍加判断就行了

欢迎光临 Fortran Coder (http://bbs.fcode.cn/)