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

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] syntaxhighlighter_viewsource syntaxhighlighter_copycode

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

kerb

你手头有maple，直接解析求出(\xi,\eta)的表达式，然后稍加判断就行了

pasuka

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

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

hipeilei

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

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

hipeilei

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

kerb

另外用迭代法求解非线性方程组应该是这样
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.

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

也就是：
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)

你检查一下是否搞混了

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的一元二次方程

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

好像是有限元等参变换，似乎你搞反了，实际上是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

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

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