坐标变换程序运行得不到结果

hipeilei

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

既然是迭代，那m n 应该有初值。

hipeilei

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

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

li913

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

fcode

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

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

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

hipeilei

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

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

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

hipeilei

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

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

kerb

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

hipeilei

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

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

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

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