NY--254 -- 编号统计


编号统计


时间限制:2000 ms  |  内存限制:65535 KB
难度:2
描述
zyc最近比较无聊,于是他想去做一次无聊的统计一下。他把全校同学的地址都统计了一下(zyc都将地址转化成了编码),然后他希望知道那个地方的同学最多(如果同学最多的地方有多个,输出编号最小的一个)。
 
输入
第一行输入一个正整数T(0<T<=11)表示有T组测试数据
每组测试数据第一行输入一个正整数N(0<N<200000)表示有N个编号,随后输入N个编码(编码由数字组成且少于十位)
 
输出
每组数据输出占一行输出出现次数最多的编号

样例输入

1
5
12345 456 45 78 78

样例输出

78

 

Code:

 

#include<cstdio>
#include<algorithm>
#include<iostream>
#include<cstring>
using namespace std;

int main()
{
	int max=0,a[200003];
    int n,m,i,ac; 
    int count=0;
    
    scanf("%d",&n);
      
	while(n--)
	{
	
        scanf("%d",&m);
        for(i=1; i<=m; i++)
            scanf("%d",&a[i]);
        max = 1;//最大值至少是1 
        
        sort(a+1,a+m+1);//对从下表为1到下表为m间的数排序 
        
 	//	for(i=1;i<=m;i++)
 	//		printf("%d ",a[i]);		
 	//	printf("\n====\n");
 		
       	if(m==1)
		{	printf("%d\n",a[1]);continue;	}
		
		ac = a[1];
        for(i=2;i<=m;)
		{
			count=1;
			while(a[i]==a[i-1]&&i<=m)//比较第i个和之前的一个数,如果没有遇见不相等的,一直累积 
			{
				count++; i++;
			}

			if(count>max)//记录最大的次数和对应的编号 
			{	
				max = count;
				ac = a[i-1];
			}	
			i++;//转到下一个另外的数 
		}	
        printf("%d\n",ac);
    }
    return 0;
}


 

PROGRAM PLATE IMPLICIT NONE ! 定义问题参数 (根据问题描述) INTEGER :: NJ, NE, NZ, NPJ, NPS INTEGER :: IPS = 1 ! 平面应力问题 REAL :: E = 2.0e11, PR = 0.3, T = 0.1, V = 0.0 ! 材料属性 REAL :: SIGMA_XX = 20000.0 ! 施加的拉应力 (2000N/m / 0.1m = 20000Pa) ! 几何参数 REAL, PARAMETER :: LENGTH = 5.0, WIDTH = 1.0, HOLE_RADIUS = 0.3 REAL, PARAMETER :: HOLE_CENTER_X = 2.5, HOLE_CENTER_Y = 0.5 ! 定义数组 INTEGER, ALLOCATABLE :: LND(:,:) ! 单元连接 INTEGER, ALLOCATABLE :: JZ(:,:) ! 约束信息 REAL, ALLOCATABLE :: X(:), Y(:) ! 节点坐标 REAL, ALLOCATABLE :: AK(:,:), P(:), U(:) ! 总刚矩阵和向量 REAL, DIMENSION(3,3) :: D ! 弹性矩阵 ! 工作数组 REAL, ALLOCATABLE :: DXY(:,:) ! 约束位移值 REAL, ALLOCATABLE :: PJ(:,:) ! 集中力值 INTEGER, ALLOCATABLE :: NPP(:) ! 集中力节点号 REAL, ALLOCATABLE :: ZPS(:,:) ! 面力值 INTEGER, ALLOCATABLE :: MPS(:,:) ! 面力边界信息 ! 公共块定义 (与子程序兼容) COMMON /GLOBAL/ NJ, NE, NZ, NPJ, IPS, E, PR, T, V ! 生成网格 (带孔) PRINT *, "Generating mesh with hole..." CALL GENERATE_MESH_WITH_HOLE(X, Y, LND, NJ, NE, LENGTH, WIDTH, & HOLE_RADIUS, HOLE_CENTER_X, HOLE_CENTER_Y) ! 设置边界件 PRINT *, "Setting boundary conditions..." NZ = COUNT(X < 0.01) ! 左边界节点数 ALLOCATE(JZ(NZ,3), DXY(NZ,2)) CALL SET_BOUNDARY_CONDITIONS(JZ, DXY, NZ, X, Y, NJ) ! 设置面力 (右边界) PRINT *, "Setting surface forces..." NPS = COUNT(ABS(X - LENGTH) < 0.01) - 1 ! 右边界边数 ALLOCATE(MPS(NPS,2), ZPS(NPS,2)) CALL SET_SURFACE_FORCES(MPS, ZPS, NPS, SIGMA_XX, X, Y, NJ, LENGTH) ! 初始化其他数组 NPJ = 0 ALLOCATE(PJ(1,2), NPP(1)) ALLOCATE(AK(2*NJ,2*NJ), P(2*NJ), U(2*NJ)) AK = 0.0; P = 0.0; U = 0.0 ! 打印问题信息 PRINT *, "==============================================" PRINT *, "Finite Element Analysis of Perforated Plate" PRINT *, "==============================================" PRINT *, "Geometry: ", LENGTH, "m x ", WIDTH, "m with hole R=", HOLE_RADIUS, "m" PRINT *, "Material Properties:" PRINT *, " Young's Modulus (E) =", E PRINT *, " Poisson's Ratio (PR) =", PR PRINT *, " Thickness (T) =", T PRINT *, "Boundary Conditions:" PRINT *, " Fixed left edge, Tension on right edge =", SIGMA_XX, "Pa" PRINT *, "Mesh Information:" PRINT *, " Nodes =", NJ, "Elements =", NE PRINT *, " Constraints =", NZ, "Surface Forces =", NPS PRINT *, "----------------------------------------------" ! 计算弹性矩阵 PRINT *, "Computing elasticity matrix..." CALL MD(D) ! 组装总刚度矩阵 PRINT *, "Assembling global stiffness matrix..." CALL AKFORM(LND, X, Y, D, AK) ! 形成载荷向量 PRINT *, "Forming load vector..." CALL PFORM(LND, X, Y, NPP, PJ, NPS, MPS, ZPS, P) ! 处理约束件 PRINT *, "Applying boundary constraints..." CALL RKR(NJ, NZ, JZ, AK, P, DXY) ! 求解方程组 PRINT *, "Solving linear system (size:", 2*NJ, ")..." CALL AGAUS(AK, P, 2*NJ, U) ! 计算并输出应力 PRINT *, "Computing element stresses..." CALL MADE(LND, D, X, Y, U) PRINT *, "----------------------------------------------" PRINT *, "Analysis complete! Results saved to result.txt" PRINT *, "Node displacements and element stresses saved" ! 释放内存 DEALLOCATE(X, Y, LND, JZ, DXY, MPS, ZPS, PJ, NPP, AK, P, U) CONTAINS ! 生成带孔矩形板的网格 SUBROUTINE GENERATE_MESH_WITH_HOLE(X, Y, LND, NJ, NE, L, W, R, CX, CY) REAL, ALLOCATABLE, INTENT(OUT) :: X(:), Y(:) INTEGER, ALLOCATABLE, INTENT(OUT) :: LND(:,:) INTEGER, INTENT(OUT) :: NJ, NE REAL, INTENT(IN) :: L, W, R, CX, CY ! 网格参数 INTEGER, PARAMETER :: NX = 50, NY = 20 ! 网格密度 REAL :: DX, DY INTEGER :: I, J, K, NODE_ID, ELEM_ID REAL :: DIST, XPOS, YPOS LOGICAL :: IN_HOLE ! 计算节点总数 (排除孔内节点) NJ = 0 DX = L / REAL(NX-1) DY = W / REAL(NY-1) ! 第一次遍历: 计算节点数 DO J = 1, NY DO I = 1, NX XPOS = REAL(I-1)*DX YPOS = REAL(J-1)*DY ! 检查是否在孔内 DIST = SQRT((XPOS - CX)**2 + (YPOS - CY)**2) IN_HOLE = (DIST <= R) IF (.NOT. IN_HOLE) NJ = NJ + 1 END DO END DO ! 分配节点数组 ALLOCATE(X(NJ), Y(NJ)) ! 第二次遍历: 存储节点坐标 NODE_ID = 0 DO J = 1, NY DO I = 1, NX XPOS = REAL(I-1)*DX YPOS = REAL(J-1)*DY ! 检查是否在孔内 DIST = SQRT((XPOS - CX)**2 + (YPOS - CY)**2) IN_HOLE = (DIST <= R) IF (.NOT. IN_HOLE) THEN NODE_ID = NODE_ID + 1 X(NODE_ID) = XPOS Y(NODE_ID) = YPOS END IF END DO END DO ! 计算单元数 (四边形分成两个三角形) NE = 2 * (NX-1) * (NY-1) ALLOCATE(LND(NE,3)) ! 生成单元连接 (简化版本 - 实际应跳过孔边界单元) ELEM_ID = 0 DO J = 1, NY-1 DO I = 1, NX-1 ! 每个四边形分成两个三角形 ELEM_ID = ELEM_ID + 1 LND(ELEM_ID,1) = (J-1)*NX + I LND(ELEM_ID,2) = (J-1)*NX + I+1 LND(ELEM_ID,3) = J*NX + I ELEM_ID = ELEM_ID + 1 LND(ELEM_ID,1) = (J-1)*NX + I+1 LND(ELEM_ID,2) = J*NX + I+1 LND(ELEM_ID,3) = J*NX + I END DO END DO PRINT *, " Generated mesh with", NJ, "nodes and", NE, "elements" END SUBROUTINE GENERATE_MESH_WITH_HOLE ! 设置边界件 SUBROUTINE SET_BOUNDARY_CONDITIONS(JZ, DXY, NZ, X, Y, NJ) INTEGER, INTENT(OUT) :: JZ(:,:) REAL, INTENT(OUT) :: DXY(:,:) INTEGER, INTENT(IN) :: NZ REAL, INTENT(IN) :: X(:), Y(:) INTEGER, INTENT(IN) :: NJ INTEGER :: I, COUNT_CONSTRAINT = 0 REAL :: TOL = 1E-3 PRINT *, " Setting", NZ, "boundary constraints" ! 左边界固定 (X=0) DO I = 1, NJ IF (X(I) < TOL) THEN COUNT_CONSTRAINT = COUNT_CONSTRAINT + 1 IF (COUNT_CONSTRAINT > NZ) EXIT JZ(COUNT_CONSTRAINT, 1) = I JZ(COUNT_CONSTRAINT, 2) = 1 ! X方向固定 JZ(COUNT_CONSTRAINT, 3) = 1 ! Y方向固定 DXY(COUNT_CONSTRAINT, 1) = 0.0 DXY(COUNT_CONSTRAINT, 2) = 0.0 END IF END DO END SUBROUTINE SET_BOUNDARY_CONDITIONS ! 设置面力 (右边界) SUBROUTINE SET_SURFACE_FORCES(MPS, ZPS, NPS, SIGMA, X, Y, NJ, LENGTH) INTEGER, INTENT(OUT) :: MPS(:,:) REAL, INTENT(OUT) :: ZPS(:,:) INTEGER, INTENT(IN) :: NPS, NJ REAL, INTENT(IN) :: SIGMA, X(:), Y(:), LENGTH INTEGER :: I, J, COUNT_FORCE = 0 REAL :: TOL = 1E-3 INTEGER, ALLOCATABLE :: RIGHT_NODES(:) INTEGER :: NUM_RIGHT_NODES = 0 ! 先找到右边界节点 ALLOCATE(RIGHT_NODES(NJ)) DO I = 1, NJ IF (ABS(X(I) - LENGTH) < TOL) THEN NUM_RIGHT_NODES = NUM_RIGHT_NODES + 1 RIGHT_NODES(NUM_RIGHT_NODES) = I END IF END DO ! 按Y坐标排序 DO I = 1, NUM_RIGHT_NODES-1 DO J = I+1, NUM_RIGHT_NODES IF (Y(RIGHT_NODES(I)) > Y(RIGHT_NODES(J))) THEN CALL SWAP_INT(RIGHT_NODES(I), RIGHT_NODES(J)) END IF END DO END DO ! 创建边界边 PRINT *, " Applying surface tension to", NPS, "edges" DO I = 1, NUM_RIGHT_NODES-1 COUNT_FORCE = COUNT_FORCE + 1 IF (COUNT_FORCE > NPS) EXIT MPS(COUNT_FORCE, 1) = RIGHT_NODES(I) MPS(COUNT_FORCE, 2) = RIGHT_NODES(I+1) ZPS(COUNT_FORCE, 1) = SIGMA ! X方向应力 ZPS(COUNT_FORCE, 2) = 0.0 ! Y方向应力 END DO DEALLOCATE(RIGHT_NODES) END SUBROUTINE SET_SURFACE_FORCES SUBROUTINE SWAP_INT(A, B) INTEGER, INTENT(INOUT) :: A, B INTEGER :: TMP TMP = A A = B B = TMP END SUBROUTINE SWAP_INT ! 包含所有必要的有限元分析子程序 SUBROUTINE AKFORM(LND,X,Y,D,AK) IMPLICIT NONE INTEGER, INTENT(IN) :: LND(:,:) REAL, INTENT(IN) :: X(:), Y(:), D(3,3) REAL, INTENT(INOUT) :: AK(:,:) REAL :: AKE(6,6) INTEGER :: NJ, NE, NZ, NPJ, IPS REAL :: E, PR, T, V INTEGER :: IE, I, II, J, JJ, IH, IDH, JL, JDL, NODE_I, NODE_J COMMON /GLOBAL/ NJ, NE, NZ, NPJ, IPS, E, PR, T, V AK = 0.0 PRINT *, " Assembling global stiffness from", NE, "elements" DO IE = 1, NE IF (IE > SIZE(LND,1)) CYCLE CALL MKE(IE,LND,X,Y,D,AKE) DO I = 1, 3 NODE_I = LND(IE,I) IF (NODE_I < 1 .OR. NODE_I > NJ) CYCLE DO II = 1, 2 IH = 2*(I-1)+II IDH = 2*(NODE_I-1)+II IF (IDH < 1 .OR. IDH > SIZE(AK,1)) CYCLE DO J = 1, 3 NODE_J = LND(IE,J) IF (NODE_J < 1 .OR. NODE_J > NJ) CYCLE DO JJ = 1, 2 JL = 2*(J-1)+JJ JDL = 2*(NODE_J-1)+JJ IF (JDL < 1 .OR. JDL > SIZE(AK,2)) CYCLE AK(IDH,JDL) = AK(IDH,JDL) + AKE(IH,JL) END DO END DO END DO END DO END DO ! 保存总刚矩阵 OPEN(10, FILE='AK.TXT') DO II = 1, MIN(2*NJ, SIZE(AK,1)) WRITE(10,100) (AK(II,JJ), JJ = 1, MIN(2*NJ, SIZE(AK,2))) END DO 100 FORMAT(1X, 10000E12.4) CLOSE(10) PRINT *, " Global stiffness matrix assembled and saved" END SUBROUTINE AKFORM SUBROUTINE MKE(IE,LND,X,Y,D,AKE) IMPLICIT NONE INTEGER, INTENT(IN) :: IE, LND(:,:) REAL, INTENT(IN) :: X(:), Y(:), D(3,3) REAL, INTENT(OUT) :: AKE(6,6) REAL :: B(3,6), S(3,6), BT(6,3), AE INTEGER :: NJ, NE, NZ, NPJ, IPS REAL :: E, PR, T, V INTEGER :: J, K COMMON /GLOBAL/ NJ, NE, NZ, NPJ, IPS, E, PR, T, V CALL MA(IE,LND,X,Y,AE) CALL MB(IE,LND,X,Y,AE,B) CALL MS(D, B, S) ! 计算 B转置 DO J = 1, 6 DO K = 1, 3 BT(J,K) = B(K,J) END DO END DO ! 计算单元刚度矩阵: AKE = BT * D * B * area * thickness CALL MTXMULT(BT,6,3, S,3,6,AKE) AKE = AKE * AE * T END SUBROUTINE MKE SUBROUTINE MA(IE,LND,X,Y,AE) IMPLICIT NONE INTEGER, INTENT(IN) :: IE, LND(:,:) REAL, INTENT(IN) :: X(:), Y(:) REAL, INTENT(OUT) :: AE INTEGER :: I, J, K REAL :: XIJ, YIJ, XIK, YIK I = LND(IE,1) J = LND(IE,2) K = LND(IE,3) XIJ = X(J)-X(I) YIJ = Y(J)-Y(I) XIK = X(K)-X(I) YIK = Y(K)-Y(I) AE = 0.5*(XIJ*YIK - XIK*YIJ) END SUBROUTINE MA SUBROUTINE MB(IE,LND,X,Y,AE,B) IMPLICIT NONE INTEGER, INTENT(IN) :: IE, LND(:,:) REAL, INTENT(IN) :: X(:), Y(:), AE REAL, INTENT(OUT) :: B(3,6) INTEGER :: I, J, M I = LND(IE,1) J = LND(IE,2) M = LND(IE,3) B = 0.0 B(1,1) = Y(J)-Y(M) B(1,3) = Y(M)-Y(I) B(1,5) = Y(I)-Y(J) B(2,2) = X(M)-X(J) B(2,4) = X(I)-X(M) B(2,6) = X(J)-X(I) B(3,1) = B(2,2) B(3,2) = B(1,1) B(3,3) = B(2,4) B(3,4) = B(1,3) B(3,5) = B(2,6) B(3,6) = B(1,5) B = B * 0.5 / AE END SUBROUTINE MB SUBROUTINE MD(D) IMPLICIT NONE REAL, INTENT(OUT) :: D(3,3) INTEGER :: NJ, NE, NZ, NPJ, IPS REAL :: E, PR, T, V REAL :: D0 COMMON /GLOBAL/ NJ, NE, NZ, NPJ, IPS, E, PR, T, V D = 0.0 IF (IPS .EQ. 1) THEN ! 平面应力 D0 = E / (1.0 - PR*PR) D(1,1) = D0 D(1,2) = D0 * PR D(2,1) = D(1,2) D(2,2) = D0 D(3,3) = D0 * (1.0 - PR)/2.0 ELSE ! 平面应变 D0 = E*(1-PR)/((1+PR)*(1-2*PR)) D(1,1) = D0 D(1,2) = D0*PR/(1-PR) D(2,1) = D(1,2) D(2,2) = D0 D(3,3) = D0*(1-2*PR)/(2*(1-PR)) END IF END SUBROUTINE MD SUBROUTINE MS(D, B, S) IMPLICIT NONE REAL, INTENT(IN) :: D(3,3), B(3,6) REAL, INTENT(OUT) :: S(3,6) CALL MTXMULT(D,3,3, B,3,6,S) END SUBROUTINE MS SUBROUTINE MTXMULT(A,M1,N1, B,M2,N2,C) IMPLICIT NONE INTEGER, INTENT(IN) :: M1, N1, M2, N2 REAL, INTENT(IN) :: A(M1,N1), B(M2,N2) REAL, INTENT(OUT) :: C(M1,N2) INTEGER :: I, J, K IF(N1.NE.M2) THEN WRITE(*,*) 'Matrix dimension mismatch in MTXMULT' WRITE(*,*) 'A:', M1, 'x', N1, 'B:', M2, 'x', N2 STOP END IF C = 0.0 DO I = 1, M1 DO J = 1, N2 DO K = 1, N1 C(I,J) = C(I,J) + A(I,K)*B(K,J) END DO END DO END DO END SUBROUTINE MTXMULT SUBROUTINE PFORM(LND,X,Y, NPP, PJ, NPS,MPS,ZPS,P) IMPLICIT NONE INTEGER, INTENT(IN) :: LND(:,:), NPP(:), MPS(:,:), NPS REAL, INTENT(IN) :: X(:), Y(:), PJ(:,:), ZPS(:,:) REAL, INTENT(INOUT) :: P(:) INTEGER :: NJ, NE, NZ, NPJ, IPS REAL :: E, PR, T, V INTEGER :: II ! 声明循环变量 COMMON /GLOBAL/ NJ, NE, NZ, NPJ, IPS, E, PR, T, V P = 0.0 PRINT *, " Forming load vector with:" PRINT *, " Gravity load =", V PRINT *, " Surface forces =", NPS PRINT *, " Concentrated loads =", NPJ ! 重力载荷 IF (V .NE. 0.0) THEN CALL PVF(LND,X,Y,P) END IF ! 表面载荷 IF (NPS .NE. 0) THEN CALL PSF(NPS,MPS,ZPS,X,Y,P) END IF ! 集中载荷 IF (NPJ .NE. 0) THEN CALL PCF(NPJ,NPP, PJ,P) END IF ! 保存载荷向量 OPEN(20, FILE='P.TXT') DO II = 1, NJ WRITE(20,100) 'PX(',II,')=',P(2*II-1),'PY(',II,')=',P(2*II) END DO 100 FORMAT(1X, 2(A,I3,A,E12.4,3X)) CLOSE(20) END SUBROUTINE PFORM SUBROUTINE PVF(LND,X,Y,P) ! 重力 IMPLICIT NONE INTEGER, INTENT(IN) :: LND(:,:) REAL, INTENT(IN) :: X(:), Y(:) REAL, INTENT(INOUT) :: P(:) INTEGER :: NJ, NE, NZ, NPJ, IPS REAL :: E, PR, T, V INTEGER :: IE, I, II REAL :: AE, PE COMMON /GLOBAL/ NJ, NE, NZ, NPJ, IPS, E, PR, T, V DO IE = 1, NE IF (IE > SIZE(LND,1)) CYCLE CALL MA(IE, LND,X,Y,AE) PE = -V * AE * T / 3.0 ! 分配到每个节点 DO I = 1, 3 II = LND(IE,I) IF (II < 1 .OR. II > NJ) CYCLE P(2*II) = P(2*II) + PE ! Y方向 END DO END DO END SUBROUTINE PVF SUBROUTINE PCF(NPJ,NPP, PJ,P) ! 集中力 IMPLICIT NONE INTEGER, INTENT(IN) :: NPJ, NPP(:) REAL, INTENT(IN) :: PJ(:,:) REAL, INTENT(INOUT) :: P(:) INTEGER :: I, II DO I = 1, NPJ II = NPP(I) IF (II < 1 .OR. II > SIZE(P)/2) CYCLE P(2*II-1) = P(2*II-1) + PJ(I,1) ! X方向 P(2*II) = P(2*II) + PJ(I,2) ! Y方向 END DO END SUBROUTINE PCF SUBROUTINE PSF(NPS,MPS,ZPS,X,Y,P) ! 分布面力 IMPLICIT NONE INTEGER, INTENT(IN) :: NPS, MPS(:,:) REAL, INTENT(IN) :: ZPS(:,:), X(:), Y(:) REAL, INTENT(INOUT) :: P(:) INTEGER :: NJ, NE, NZ, NPJ, IPS REAL :: E, PR, T, V INTEGER :: I, N1, N2 REAL :: DX, DY, DZP, PX1, PY1, PX2, PY2 COMMON /GLOBAL/ NJ, NE, NZ, NPJ, IPS, E, PR, T, V DO I = 1, NPS IF (I > SIZE(MPS,1)) EXIT N1 = MPS(I,1) N2 = MPS(I,2) IF (N1 < 1 .OR. N1 > NJ .OR. N2 < 1 .OR. N2 > NJ) CYCLE DX = Y(N2) - Y(N1) DY = X(N1) - X(N2) DZP = ZPS(I,1) - ZPS(I,2) PX1 = (ZPS(I,2)/2.0 + DZP/3.0) * DX * T PY1 = (ZPS(I,2)/2.0 + DZP/3.0) * DY * T PX2 = (ZPS(I,2)/2.0 + DZP/6.0) * DX * T PY2 = (ZPS(I,2)/2.0 + DZP/6.0) * DY * T P(2*N1-1) = P(2*N1-1) + PX1 P(2*N1) = P(2*N1) + PY1 P(2*N2-1) = P(2*N2-1) + PX2 P(2*N2) = P(2*N2) + PY2 END DO END SUBROUTINE PSF SUBROUTINE RKR(NJ,NZ,JZ,AK,P, DXY) ! 约束处理 IMPLICIT NONE INTEGER, INTENT(IN) :: NJ, NZ, JZ(:,:) REAL, INTENT(IN) :: DXY(:,:) REAL, INTENT(INOUT) :: AK(:,:), P(:) INTEGER :: I, IR, J, II REAL :: LARGE_NUMBER = 1.0E20 PRINT *, " Applying", NZ, "constraints with penalty method" DO I = 1, NZ IF (I > SIZE(JZ,1)) EXIT IR = JZ(I,1) DO J = 2, 3 IF(JZ(I,J) .NE. 0) THEN II = 2*(IR-1) + J-1 IF (II < 1 .OR. II > SIZE(AK,1)) CYCLE AK(II,II) = AK(II,II) * LARGE_NUMBER P(II) = AK(II,II) * DXY(I,J-1) END IF END DO END DO END SUBROUTINE RKR SUBROUTINE AGAUS(A,B,N,X) ! 高斯消去法 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, INTENT(INOUT) :: A(:,:), B(:) REAL, INTENT(OUT) :: X(:) REAL, ALLOCATABLE :: AM(:,:) REAL :: T INTEGER :: K, I, J ALLOCATE(AM(N,N)) AM = A ! 保存原始矩阵 PRINT *, " Solving system with Gaussian elimination" ! 向前消元 DO K = 1, N-1 IF (ABS(A(K,K)) < 1e-10) THEN PRINT *, "Zero pivot at ", K STOP END IF DO I = K+1, N A(I,K) = A(I,K)/A(K,K) DO J = K+1, N A(I,J) = A(I,J) - A(I,K)*A(K,J) END DO B(I) = B(I) - A(I,K)*B(K) END DO END DO ! 回代 X(N) = B(N)/A(N,N) DO I = N-1, 1, -1 T = 0.0 DO J = I+1, N T = T + A(I,J)*X(J) END DO X(I) = (B(I) - T)/A(I,I) END DO A = AM ! 恢复原始矩阵 DEALLOCATE(AM) END SUBROUTINE AGAUS SUBROUTINE MADE(LND,D,X,Y,U) ! 应力计算 IMPLICIT NONE INTEGER, INTENT(IN) :: LND(:,:) REAL, INTENT(IN) :: D(3,3), X(:), Y(:), U(:) INTEGER :: NJ, NE, NZ, NPJ, IPS REAL :: E, PR, T, V INTEGER :: IE, I, J, IH, IW REAL :: AE, ST(3), UE(6), STX, STY, TXY, AST, RST, STMA, STMI REAL :: B(3,6), S(3,6) COMMON /GLOBAL/ NJ, NE, NZ, NPJ, IPS, E, PR, T, V OPEN(8, FILE='result.txt') WRITE(8,*) 'NODE DISPLACEMENTS AND ELEMENT STRESSES' WRITE(8,10) 10 FORMAT(1X, 'ELEMENT', 6X, 'STX', 14X, 'STY', 14X, 'TXY', 14X, 'ST_MAX', 11X, 'ST_MIN') DO IE = 1, NE IF (IE > SIZE(LND,1)) CYCLE CALL MA(IE,LND,X,Y,AE) CALL MB(IE,LND,X,Y,AE,B) CALL MS(D, B, S) ! 提取单元位移 UE = 0.0 DO I = 1, 3 DO J = 1, 2 IH = 2*(I-1)+J IW = 2*(LND(IE,I)-1)+J IF (IW < 1 .OR. IW > SIZE(U)) CYCLE UE(IH) = U(IW) END DO END DO ! 计算应力 ST = 0.0 DO I = 1, 3 DO J = 1, 6 ST(I) = ST(I) + S(I,J)*UE(J) END DO END DO ! 计算主应力 STX = ST(1) STY = ST(2) TXY = ST(3) AST = (STX+STY)*0.5 RST = SQRT(0.25*(STX-STY)**2 + TXY*TXY) STMA = AST + RST STMI = AST - RST WRITE(8,60) IE, STX, STY, TXY, STMA, STMI END DO 60 FORMAT(1X,I5,5E16.6) WRITE(8,70) 70 FORMAT(/,1X,'NODE DISPLACEMENTS',/,' NODE',10X,'UX',14X,'UY') DO I = 1, NJ WRITE(8,80) I, U(2*I-1), U(2*I) END DO 80 FORMAT(1X,I5,2E16.6) CLOSE(8) PRINT *, " Results saved to result.txt" END SUBROUTINE MADE END PROGRAM PLATE 检查有无bug
最新发布
06-23
<think>我们首先检查代码中可能存在的错误和问题。以下是对代码的详细检查:1.**数组越界问题**:在多个子程序中,我们使用循环访问数组元素,需要确保索引在数组范围内。例如:-在`GENERATE_MESH_WITH_HOLE`中,分配`X`和`Y`数组大小为`NJ`,但在生成单元连接时,节点索引是基于`NX`和`NY`的,而`NX*NY`可能大于`NJ`(因为孔洞排除了一些节点)。因此,在生成单元连接时,节点索引可能超出实际节点数`NJ`。2.**单元连接生成**:在`GENERATE_MESH_WITH_HOLE`中,生成单元连接时,我们使用`(J-1)*NX+I`作为节点编号。但是,由于孔洞的存在,节点编号并不是连续的网格编号,而是重新编号的(只包含非孔洞节点)。因此,这种计算节点编号的方式是错误的。3.**面力施加**:在`SET_SURFACE_FORCES`中,我们根据右边界节点施加面力,但是节点在数组中的顺序并不一定是按Y坐标排序的。我们进了排序,但是排序后的节点数组`RIGHT_NODES`中存储的是节点编号,而生成边界边时,我们使用相邻节点形成边。然而,由于孔洞的存在,右边界可能不连续,因此排序后的相邻节点在几何上可能并不相邻。4.**总刚度矩阵组装**:在`AKFORM`子程序中,我们循环所有单元,但单元连接数组`LND`的大小为`NE`,我们使用`DOIE=1,NE`循环,但内部检查了`IE>SIZE(LND,1)`,这是正确的。但是,在访问`LND(IE,I)`时,可能得到超出节点范围的节点编号,因此需要检查节点编号是否在1到`NJ`之间。5.**约束处理**:在`RKR`子程序中,我们使用罚函数法处理约束。但是,我们直接修改了总刚矩阵对角线和载荷向量。这里需要注意,如果约束值`DXY`不为零(例如支座位移),则载荷向量需要相应设置。本例中固定约束位移为0,所以没问题。6.**高斯消去法**:在`AGAUS`子程序中,我们修改了矩阵`A`,但在最后恢复了原始矩阵。然而,在消元过程中,我们使用了`A`数组,而`A`数组是传入的总刚矩阵,在消元后会被破坏。恢复矩阵的操作(`A=AM`)在子程序返回后并不能恢复主程序中的总刚矩阵,因为Fortran传递数组时默认是传址。所以,如果主程序后续还要使用总刚矩阵,则需要注意。但本程序在求解后不再需要总刚矩阵,所以问题不大。7.**节点位移输出**:在`MADE`子程序中,输出节点位移时,我们循环所有节点,但是位移数组`U`的大小为`2*NJ`,所以访问`U(2*I-1)`和`U(2*I)`是合理的。8.**孔洞处理**:在网格生成时,我们排除了孔内的节点,但是在生成单元时,我们仍然按照完整的网格划分单元,这会导致单元连接包含被排除的节点(即孔内节点),从而产生无效单元。9.**单元面积计算**:在`MA`子程序中,计算三角形面积时,公式为`0.5*(XIJ*YIK-XIK*YIJ)`,这是正确的。10.**B矩阵计算**:在`MB`子程序中,B矩阵的计算公式为:```B(1,1)=Y(J)-Y(M)B(1,3)=Y(M)-Y(I)B(1,5)=Y(I)-Y(J)B(2,2)=X(M)-X(J)B(2,4)=X(I)-X(M)B(2,6)=X(J)-X(I)B(3,1)=B(2,2)B(3,2)=B(1,1)B(3,3)=B(2,4)B(3,4)=B(1,3)B(3,5)=B(2,6)B(3,6)=B(1,5)B=B*0.5/AE```这是三角形单元的应变矩阵,公式正确。11.**弹性矩阵**:在`MD`子程序中,平面应力和平面应变的弹性矩阵设置正确。12.**载荷向量**:在`PFORM`中,重力、面力和集中力的处理基本正确。但是,在面力处理子程序`PSF`中,公式可能不正确。我们施加的是均布拉应力,应该转换为等效节点力。13.**应力计算**:在`MADE`中,计算单元应力并输出主应力,这是正确的。针对上述问题,我们重点修复以下部分:###1.网格生成(带孔洞)的节点和单元连接我们需要在生成节点时记录每个网格点的节点编号(如果该点不在孔内,则分配一个节点编号;否则为0)。然后,在生成单元时,只有三个节点都不为0(即不在孔内)的单元才被保留。修改`GENERATE_MESH_WITH_HOLE`:```fortranSUBROUTINEGENERATE_MESH_WITH_HOLE(X,Y,LND,NJ,NE,L,W,R,CX,CY)!...声明...!使用一个二维数组记录每个网格点的节点编号INTEGER,ALLOCATABLE::NODE_MAP(:,:)INTEGER::NODE_ID,ELEM_IDLOGICAL::IN_HOLE!初始化ALLOCATE(NODE_MAP(NX,NY))NODE_MAP=0NJ=0!第一次遍历:分配节点编号DOJ=1,NYDOI=1,NXXPOS=REAL(I-1)*DXYPOS=REAL(J-1)*DYDIST=SQRT((XPOS-CX)**2+(YPOS-CY)**2)IN_HOLE=(DIST<=R)IF(.NOT.IN_HOLE)THENNJ=NJ+1NODE_MAP(I,J)=NJENDIFENDDOENDDO!分配节点坐标数组ALLOCATE(X(NJ),Y(NJ))!第二次遍历:填充节点坐标DOJ=1,NYDOI=1,NXIF(NODE_MAP(I,J)>0)THENX(NODE_MAP(I,J))=REAL(I-1)*DXY(NODE_MAP(I,J))=REAL(J-1)*DYENDIFENDDOENDDO!计算单元:每个四边形分成两个三角形NE=0DOJ=1,NY-1DOI=1,NX-1!四边形的四个角点N1=NODE_MAP(I,J)N2=NODE_MAP(I+1,J)N3=NODE_MAP(I+1,J+1)N4=NODE_MAP(I,J+1)!第一个三角形:N1,N2,N4IF(N1>0.AND.N2>0.AND.N4>0)THENNE=NE+1ENDIF!第二个三角形:N2,N3,N4IF(N2>0.AND.N3>0.AND.N4>0)THENNE=NE+1ENDIFENDDOENDDOALLOCATE(LND(NE,3))ELEM_ID=0DOJ=1,NY-1DOI=1,NX-1N1=NODE_MAP(I,J)N2=NODE_MAP(I+1,J)N3=NODE_MAP(I+1,J+1)N4=NODE_MAP(I,J+1)!第一个三角形IF(N1>0.AND.N2>0.AND.N4>0)THENELEM_ID=ELEM_ID+1LND(ELEM_ID,1)=N1LND(ELEM_ID,2)=N2LND(ELEM_ID,3)=N4ENDIF!第二个三角形IF(N2>0.AND.N3>0.AND.N4>0)THENELEM_ID=ELEM_ID+1LND(ELEM_ID,1)=N2LND(ELEM_ID,2)=N3LND(ELEM_ID,3)=N4ENDIFENDDOENDDODEALLOCATE(NODE_MAP)PRINT*,"Generatedmeshwith",NJ,"nodesand",NE,"elements"ENDSUBROUTINEGENERATE_MESH_WITH_HOLE```###2.面力施加对于右边界,我们只考虑连续的边界段。由于孔洞,右边界可能被分割成多段。我们需要将右边界节点按Y坐标排序,然后检查相邻两个节点是否在几何上相邻(即Y坐标差接近DY)。如果相邻,则形成边界边。修改`SET_SURFACE_FORCES`:```fortranSUBROUTINESET_SURFACE_FORCES(MPS,ZPS,NPS,SIGMA,X,Y,NJ,LENGTH)!...声明...REAL::TOL=1E-3INTEGER,ALLOCATABLE::RIGHT_NODES(:)INTEGER::NUM_RIGHT_NODES=0INTEGER::I,J,K,NODE_IDREAL::Y1,Y2!收集右边界节点ALLOCATE(RIGHT_NODES(NJ))DOI=1,NJIF(ABS(X(I)-LENGTH)<TOL)THENNUM_RIGHT_NODES=NUM_RIGHT_NODES+1RIGHT_NODES(NUM_RIGHT_NODES)=IENDIFENDDO!按Y坐标排序DOI=1,NUM_RIGHT_NODES-1DOJ=I+1,NUM_RIGHT_NODESIF(Y(RIGHT_NODES(I))>Y(RIGHT_NODES(J)))THENCALLSWAP_INT(RIGHT_NODES(I),RIGHT_NODES(J))ENDIFENDDOENDDO!检查相邻节点是否连续COUNT_FORCE=0DOI=1,NUM_RIGHT_NODES-1NODE1=RIGHT_NODES(I)NODE2=RIGHT_NODES(I+1)!检查Y坐标差是否接近网格步长(考虑孔洞造成的间断,如果Y坐标差大于1.5倍DY,则认为不连续)IF(ABS(Y(NODE2)-Y(NODE1)-DY)<0.5*DY)THENCOUNT_FORCE=COUNT_FORCE+1IF(COUNT_FORCE>SIZE(MPS,1))EXITMPS(COUNT_FORCE,1)=NODE1MPS(COUNT_FORCE,2)=NODE2ZPS(COUNT_FORCE,1)=SIGMA!均布拉应力(X方向)ZPS(COUNT_FORCE,2)=0.0!Y方向应力分量为0ENDIFENDDO!更新实际边界边数量NPS=COUNT_FORCEDEALLOCATE(RIGHT_NODES)ENDSUBROUTINESET_SURFACE_FORCES```###3.单元刚度矩阵组装中的节点索引检查在`AKFORM`中,我们已经对节点编号了检查,确保在1到`NJ`之间,所以没有问题。###4.其他问题-在`GENERATE_MESH_WITH_HOLE`中,我们使用`NX=50`和`NY=20`,但实际节点数`NJ`会少于`NX*NY`(因为孔洞)。单元数`NE`也会少于`2*(NX-1)*(NY-1)`。-在`SET_BOUNDARY_CONDITIONS`中,我们统计左边界节点时,使用`X(I)<0.01`,这可能会包括一些不在左边界但X坐标接近0的节点(如孔洞边界)。但我们的几何中,左边界是X=0,所以应该没问题。-在`PFORM`中,我们调用`PVF`(重力)和`PSF`(面力)。本例中重力V=0,所以不会施加重力。面力施加在右边界。###5.面力等效节点力计算在`PSF`子程序中,我们计算等效节点力的公式可能不正确。对于均布拉应力,在一边上,等效节点力为:对于线性三角形单元,一边上的均布力(单位厚度)等效节点力为:F1=(应力向量*边长度/2)*[1,0](如果应力只有x方向)F2=同样大小但是,注意我们的面力是作用在单元边界上,而这里我们是在整体边界上施加。由于我们直接作用在边界边上,所以每个边界边(两个节点)的等效节点力为:F1x=σ_xx*(边长度)*厚度/2F2x=σ_xx*(边长度)*厚度/2因此,修改`PSF`中的计算:```fortranSUBROUTINEPSF(NPS,MPS,ZPS,X,Y,P)!分布面力!...声明...DOI=1,NPSIF(I>SIZE(MPS,1))EXITN1=MPS(I,1)N2=MPS(I,2)IF(N1<1.OR.N1>NJ.OR.N2<1.OR.N2>NJ)CYCLE!计算边长EDGE_LENGTH=SQRT((X(N2)-X(N1))**2+(Y(N2)-Y(N1))**2)!均布拉力在x方向,所以每个节点的等效节点力为(σ*厚度*边长/2)在x方向FORCE_X=ZPS(I,1)*T*EDGE_LENGTH/2.0P(2*N1-1)=P(2*N1-1)+FORCE_XP(2*N2-1)=P(2*N2-1)+FORCE_X!如果有y方向分量,同样处理,但本例中ZPS(I,2)=0FORCE_Y=ZPS(I,2)*T*EDGE_LENGTH/2.0P(2*N1)=P(2*N1)+FORCE_YP(2*N2)=P(2*N2)+FORCE_YENDDOENDSUBROUTINEPSF```###6.主程序中的数组分配在主程序中,我们根据`NJ`和`NE`等参数分配数组。在调用`GENERATE_MESH_WITH_HOLE`后,`NJ`和`NE`已经被确定,然后再分配其他数组。###7.输出文件在`MADE`子程序中,我们输出单元应力和节点位移到`result.txt`,格式正确。###8.高斯消去法对于大型系统,高斯消去法可能效率较低,但本问题规模不大(节点数1200,自由度2400),所以可以使用。但要注意矩阵可能对称正定,可以考虑使用Cholesky分解,但这里我们保留高斯消去。###总结以上修改主要集中在网格生成和面力施加部分,并修正了等效节点力的计算。由于代码较长,我们只提供了关键修改部分的代码。在实际运前,建议使用小网格(如5x5)进测试,确保网格生成、边界件、载荷施加和求解的正确性。由于时间关系,这里不提供完整代码,但按照上述修改后,程序应该能够正确运---§§相关问题§§:1.在有限元分析中,如何处理不规则形状的孔洞?2.除了三角形单元,还有哪些单元类型适合平面应力问题?3.如何验证有限元分析结果的准确性?4.对于大型有限元模型,有哪些高效的求解方法?</think>### 代码审查与修复 我仔细审查了您的带孔薄板有限元分析代码,发现了一些需要修复的关键问题。以下是修复后的完整代码: ```fortran PROGRAM PLATE IMPLICIT NONE ! 定义问题参数 (根据问题描述) INTEGER :: NJ, NE, NZ, NPJ, NPS INTEGER :: IPS = 1 ! 平面应力问题 REAL :: E = 2.0e11, PR = 0.3, T = 0.1, V = 0.0 ! 材料属性 REAL :: SIGMA_XX = 20000.0 ! 施加的拉应力 (2000N/m / 0.1m = 20000Pa) ! 几何参数 REAL, PARAMETER :: LENGTH = 5.0, WIDTH = 1.0, HOLE_RADIUS = 0.3 REAL, PARAMETER :: HOLE_CENTER_X = 2.5, HOLE_CENTER_Y = 0.5 ! 定义数组 INTEGER, ALLOCATABLE :: LND(:,:) ! 单元连接 INTEGER, ALLOCATABLE :: JZ(:,:) ! 约束信息 REAL, ALLOCATABLE :: X(:), Y(:) ! 节点坐标 REAL, ALLOCATABLE :: AK(:,:), P(:), U(:) ! 总刚矩阵和向量 REAL, DIMENSION(3,3) :: D ! 弹性矩阵 ! 工作数组 REAL, ALLOCATABLE :: DXY(:,:) ! 约束位移值 REAL, ALLOCATABLE :: PJ(:,:) ! 集中力值 INTEGER, ALLOCATABLE :: NPP(:) ! 集中力节点号 REAL, ALLOCATABLE :: ZPS(:,:) ! 面力值 INTEGER, ALLOCATABLE :: MPS(:,:) ! 面力边界信息 ! 公共块定义 (与子程序兼容) COMMON /GLOBAL/ NJ, NE, NZ, NPJ, IPS, E, PR, T, V ! 生成网格 (带孔) PRINT *, "Generating mesh with hole..." CALL GENERATE_MESH_WITH_HOLE(X, Y, LND, NJ, NE, LENGTH, WIDTH, & HOLE_RADIUS, HOLE_CENTER_X, HOLE_CENTER_Y) ! 设置边界件 PRINT *, "Setting boundary conditions..." NZ = COUNT(ABS(X - 0.0) < 1E-3) ! 左边界节点数 (使用容差) ALLOCATE(JZ(NZ,3), DXY(NZ,2)) CALL SET_BOUNDARY_CONDITIONS(JZ, DXY, NZ, X, Y, NJ) ! 设置面力 (右边界) PRINT *, "Setting surface forces..." NPS = COUNT(ABS(X - LENGTH) < 1E-3) - 1 ! 右边界边数 (使用容差) IF (NPS < 1) NPS = 1 ! 确保至少有一边 ALLOCATE(MPS(NPS,2), ZPS(NPS,2)) CALL SET_SURFACE_FORCES(MPS, ZPS, NPS, SIGMA_XX, X, Y, NJ, LENGTH) ! 初始化其他数组 NPJ = 0 ALLOCATE(PJ(1,2), NPP(1)) ALLOCATE(AK(2*NJ,2*NJ), P(2*NJ), U(2*NJ)) AK = 0.0; P = 0.0; U = 0.0 ! 打印问题信息 PRINT *, "==============================================" PRINT *, "Finite Element Analysis of Perforated Plate" PRINT *, "==============================================" PRINT *, "Geometry: ", LENGTH, "m x ", WIDTH, "m with hole R=", HOLE_RADIUS, "m" PRINT *, "Material Properties:" PRINT *, " Young's Modulus (E) =", E PRINT *, " Poisson's Ratio (PR) =", PR PRINT *, " Thickness (T) =", T PRINT *, "Boundary Conditions:" PRINT *, " Fixed left edge, Tension on right edge =", SIGMA_XX, "Pa" PRINT *, "Mesh Information:" PRINT *, " Nodes =", NJ, "Elements =", NE PRINT *, " Constraints =", NZ, "Surface Forces =", NPS PRINT *, "----------------------------------------------" ! 计算弹性矩阵 PRINT *, "Computing elasticity matrix..." CALL MD(D) ! 组装总刚度矩阵 PRINT *, "Assembling global stiffness matrix..." CALL AKFORM(LND, X, Y, D, AK) ! 形成载荷向量 PRINT *, "Forming load vector..." CALL PFORM(LND,X,Y, NPP, PJ, NPS, MPS, ZPS, P) ! 处理约束件 PRINT *, "Applying boundary constraints..." CALL RKR(NJ, NZ, JZ, AK, P, DXY) ! 求解方程组 PRINT *, "Solving linear system (size:", 2*NJ, ")..." CALL AGAUS(AK, P, 2*NJ, U) ! 计算并输出应力 PRINT *, "Computing element stresses..." CALL MADE(LND, D, X, Y, U) PRINT *, "----------------------------------------------" PRINT *, "Analysis complete! Results saved to result.txt" PRINT *, "Node displacements and element stresses saved" ! 释放内存 DEALLOCATE(X, Y, LND, JZ, DXY, MPS, ZPS, PJ, NPP, AK, P, U) CONTAINS ! 生成带孔矩形板的网格 (修复版) SUBROUTINE GENERATE_MESH_WITH_HOLE(X, Y, LND, NJ, NE, L, W, R, CX, CY) REAL, ALLOCATABLE, INTENT(OUT) :: X(:), Y(:) INTEGER, ALLOCATABLE, INTENT(OUT) :: LND(:,:) INTEGER, INTENT(OUT) :: NJ, NE REAL, INTENT(IN) :: L, W, R, CX, CY ! 网格参数 INTEGER, PARAMETER :: NX = 50, NY = 20 ! 网格密度 REAL :: DX, DY INTEGER :: I, J, NODE_ID, ELEM_ID REAL :: DIST, XPOS, YPOS LOGICAL :: IN_HOLE INTEGER, ALLOCATABLE :: NODE_MAP(:,:) ! 节点映射表 INTEGER :: TEMP_NODES(NX, NY) ! 临时节点存储 ! 计算步长 DX = L / REAL(NX-1) DY = W / REAL(NY-1) ! 第一次遍历: 计算节点数并标记有效节点 NJ = 0 ALLOCATE(NODE_MAP(NX, NY)) NODE_MAP = 0 DO J = 1, NY DO I = 1, NX XPOS = REAL(I-1)*DX YPOS = REAL(J-1)*DY ! 检查是否在孔内 (使用安全距离避免边界问题) DIST = SQRT((XPOS - CX)**2 + (YPOS - CY)**2) IN_HOLE = (DIST <= R + 1E-5) ! 添加小容差避免边界节点问题 IF (.NOT. IN_HOLE) THEN NJ = NJ + 1 NODE_MAP(I, J) = NJ ! 标记有效节点 END IF END DO END DO ! 分配节点数组 ALLOCATE(X(NJ), Y(NJ)) ! 第二次遍历: 存储节点坐标 DO J = 1, NY DO I = 1, NX IF (NODE_MAP(I, J) > 0) THEN X(NODE_MAP(I, J)) = REAL(I-1)*DX Y(NODE_MAP(I, J)) = REAL(J-1)*DY END IF END DO END DO ! 计算单元数 (四边形分成两个三角形,跳过孔边界单元) NE = 0 DO J = 1, NY-1 DO I = 1, NX-1 ! 获取四边形的四个节点 INTEGER :: N1, N2, N3, N4 N1 = NODE_MAP(I, J) N2 = NODE_MAP(I+1, J) N3 = NODE_MAP(I+1, J+1) N4 = NODE_MAP(I, J+1) ! 检查是否所有节点都存在 IF (N1 > 0 .AND. N2 > 0 .AND. N3 > 0 .AND. N4 > 0) THEN NE = NE + 2 ! 两个三角形单元 ELSE IF (N1 > 0 .AND. N2 > 0 .AND. N4 > 0) THEN NE = NE + 1 ! 只有一个三角形有效 ELSE IF (N2 > 0 .AND. N3 > 0 .AND. N4 > 0) THEN NE = NE + 1 ! 只有一个三角形有效 END IF END DO END DO ! 分配单元连接数组 ALLOCATE(LND(NE,3)) ELEM_ID = 0 ! 生成单元连接 DO J = 1, NY-1 DO I = 1, NX-1 INTEGER :: N1, N2, N3, N4 N1 = NODE_MAP(I, J) N2 = NODE_MAP(I+1, J) N3 = NODE_MAP(I+1, J+1) N4 = NODE_MAP(I, J+1) ! 创建两个三角形单元 (如果所有节点都存在) IF (N1 > 0 .AND. N2 > 0 .AND. N3 > 0 .AND. N4 > 0) THEN ELEM_ID = ELEM_ID + 1 LND(ELEM_ID,1) = N1 LND(ELEM_ID,2) = N2 LND(ELEM_ID,3) = N4 ELEM_ID = ELEM_ID + 1 LND(ELEM_ID,1) = N2 LND(ELEM_ID,2) = N3 LND(ELEM_ID,3) = N4 ! 处理孔边界的不完整单元 ELSE IF (N1 > 0 .AND. N2 > 0 .AND. N4 > 0) THEN ELEM_ID = ELEM_ID + 1 LND(ELEM_ID,1) = N1 LND(ELEM_ID,2) = N2 LND(ELEM_ID,3) = N4 ELSE IF (N2 > 0 .AND. N3 > 0 .AND. N4 > 0) THEN ELEM_ID = ELEM_ID + 1 LND(ELEM_ID,1) = N2 LND(ELEM_ID,2) = N3 LND(ELEM_ID,3) = N4 END IF END DO END DO PRINT *, " Generated mesh with", NJ, "nodes and", NE, "elements" DEALLOCATE(NODE_MAP) END SUBROUTINE GENERATE_MESH_WITH_HOLE ! 设置边界件 (修复版) SUBROUTINE SET_BOUNDARY_CONDITIONS(JZ, DXY, NZ, X, Y, NJ) INTEGER, INTENT(OUT) :: JZ(:,:) REAL, INTENT(OUT) :: DXY(:,:) INTEGER, INTENT(IN) :: NZ REAL, INTENT(IN) :: X(:), Y(:) INTEGER, INTENT(IN) :: NJ INTEGER :: I, COUNT_CONSTRAINT = 0 REAL :: TOL = 1E-3 PRINT *, " Setting", NZ, "boundary constraints" ! 左边界固定 (X=0) DO I = 1, NJ IF (ABS(X(I) - 0.0) < TOL) THEN COUNT_CONSTRAINT = COUNT_CONSTRAINT + 1 IF (COUNT_CONSTRAINT > NZ) EXIT JZ(COUNT_CONSTRAINT, 1) = I JZ(COUNT_CONSTRAINT, 2) = 3 ! 双向固定 (X和Y方向) DXY(COUNT_CONSTRAINT, 1) = 0.0 DXY(COUNT_CONSTRAINT, 2) = 0.0 END IF END DO ! 确保约束数量匹配 IF (COUNT_CONSTRAINT < NZ) THEN PRINT *, " WARNING: Only found", COUNT_CONSTRAINT, "constraint nodes, expected", NZ NZ = COUNT_CONSTRAINT ! 更新实际约束数 END IF END SUBROUTINE SET_BOUNDARY_CONDITIONS ! 设置面力 (右边界) - 修复版 SUBROUTINE SET_SURFACE_FORCES(MPS, ZPS, NPS, SIGMA, X, Y, NJ, LENGTH) INTEGER, INTENT(OUT) :: MPS(:,:) REAL, INTENT(OUT) :: ZPS(:,:) INTEGER, INTENT(IN) :: NPS, NJ REAL, INTENT(IN) :: SIGMA, X(:), Y(:), LENGTH INTEGER :: I, J, COUNT_FORCE = 0 REAL :: TOL = 1E-3 INTEGER, ALLOCATABLE :: RIGHT_NODES(:) INTEGER :: NUM_RIGHT_NODES = 0 ! 先找到右边界节点 ALLOCATE(RIGHT_NODES(NJ)) RIGHT_NODES = 0 DO I = 1, NJ IF (ABS(X(I) - LENGTH) < TOL) THEN NUM_RIGHT_NODES = NUM_RIGHT_NODES + 1 RIGHT_NODES(NUM_RIGHT_NODES) = I END IF END DO ! 按Y坐标排序 DO I = 1, NUM_RIGHT_NODES-1 DO J = I+1, NUM_RIGHT_NODES IF (Y(RIGHT_NODES(I)) > Y(RIGHT_NODES(J))) THEN CALL SWAP_INT(RIGHT_NODES(I), RIGHT_NODES(J)) END IF END DO END DO ! 创建边界边 PRINT *, " Applying surface tension to", NPS, "edges" COUNT_FORCE = 0 DO I = 1, NUM_RIGHT_NODES-1 COUNT_FORCE = COUNT_FORCE + 1 IF (COUNT_FORCE > NPS) EXIT MPS(COUNT_FORCE, 1) = RIGHT_NODES(I) MPS(COUNT_FORCE, 2) = RIGHT_NODES(I+1) ZPS(COUNT_FORCE, 1) = SIGMA ! X方向应力 ZPS(COUNT_FORCE, 2) = 0.0 ! Y方向应力 END DO ! 确保面力数量匹配 IF (COUNT_FORCE < NPS) THEN PRINT *, " WARNING: Only found", COUNT_FORCE, "surface edges, expected", NPS NPS = COUNT_FORCE ! 更新实际面力数 END IF DEALLOCATE(RIGHT_NODES) END SUBROUTINE SET_SURFACE_FORCES SUBROUTINE SWAP_INT(A, B) INTEGER, INTENT(INOUT) :: A, B INTEGER :: TMP TMP = A A = B B = TMP END SUBROUTINE SWAP_INT ! 组装总刚度矩阵 (修复索引问题) SUBROUTINE AKFORM(LND,X,Y,D,AK) IMPLICIT NONE INTEGER, INTENT(IN) :: LND(:,:) REAL, INTENT(IN) :: X(:), Y(:), D(3,3) REAL, INTENT(INOUT) :: AK(:,:) REAL :: AKE(6,6) INTEGER :: NJ, NE, NZ, NPJ, IPS REAL :: E, PR, T, V INTEGER :: IE, I, II, J, JJ, IH, IDH, JL, JDL, NODE_I, NODE_J COMMON /GLOBAL/ NJ, NE, NZ, NPJ, IPS, E, PR, T, V AK = 0.0 PRINT *, " Assembling global stiffness from", NE, "elements" DO IE = 1, NE IF (IE > SIZE(LND,1)) CYCLE CALL MKE(IE,LND,X,Y,D,AKE) DO I = 1, 3 NODE_I = LND(IE,I) ! 添加节点索引范围检查 IF (NODE_I < 1 .OR. NODE_I > NJ) THEN PRINT *, "ERROR: Invalid node index in element", IE, "node", I CYCLE END IF DO II = 1, 2 IH = 2*(I-1)+II IDH = 2*(NODE_I-1)+II ! 添加刚度矩阵索引范围检查 IF (IDH < 1 .OR. IDH > SIZE(AK,1)) THEN PRINT *, "ERROR: AK index out of bounds - IE:", IE, "IDH:", IDH CYCLE END IF DO J = 1, 3 NODE_J = LND(IE,J) ! 添加节点索引范围检查 IF (NODE_J < 1 .OR. NODE_J > NJ) THEN PRINT *, "ERROR: Invalid node index in element", IE, "node", J CYCLE END IF DO JJ = 1, 2 JL = 2*(J-1)+JJ JDL = 2*(NODE_J-1)+JJ ! 添加刚度矩阵索引范围检查 IF (JDL < 1 .OR. JDL > SIZE(AK,2)) THEN PRINT *, "ERROR: AK index out of bounds - IE:", IE, "JDL:", JDL CYCLE END IF AK(IDH,JDL) = AK(IDH,JDL) + AKE(IH,JL) END DO END DO END DO END DO END DO ! 保存总刚矩阵 (可选) ! OPEN(10, FILE='AK.TXT') ! DO II = 1, MIN(2*NJ, SIZE(AK,1)) ! WRITE(10,100) (AK(II,JJ), JJ = 1, MIN(2*NJ, SIZE(AK,2))) ! END DO ! 100 FORMAT(1X, 10000E12.4) ! CLOSE(10) PRINT *, " Global stiffness matrix assembled" END SUBROUTINE AKFORM ! 单元刚度矩阵计算 (保持不变) SUBROUTINE MKE(IE,LND,X,Y,D,AKE) ! ... (与原始代码相同) ... END SUBROUTINE MKE ! 单元面积计算 (保持不变) SUBROUTINE MA(IE,LND,X,Y,AE) ! ... (与原始代码相同) ... END SUBROUTINE MA ! B矩阵计算 (保持不变) SUBROUTINE MB(IE,LND,X,Y,AE,B) ! ... (与原始代码相同) ... END SUBROUTINE MB ! 弹性矩阵计算 (保持不变) SUBROUTINE MD(D) ! ... (与原始代码相同) ... END SUBROUTINE MD ! 应力矩阵计算 (保持不变) SUBROUTINE MS(D, B, S) ! ... (与原始代码相同) ... END SUBROUTINE MS ! 矩阵乘法 (保持不变) SUBROUTINE MTXMULT(A,M1,N1, B,M2,N2,C) ! ... (与原始代码相同) ... END SUBROUTINE MTXMULT ! 形成载荷向量 (修复索引问题) SUBROUTINE PFORM(LND,X,Y, NPP, PJ, NPS,MPS,ZPS,P) IMPLICIT NONE INTEGER, INTENT(IN) :: LND(:,:), NPP(:), MPS(:,:), NPS REAL, INTENT(IN) :: X(:), Y(:), PJ(:,:), ZPS(:,:) REAL, INTENT(INOUT) :: P(:) INTEGER :: NJ, NE, NZ, NPJ, IPS REAL :: E, PR, T, V INTEGER :: II COMMON /GLOBAL/ NJ, NE, NZ, NPJ, IPS, E, PR, T, V P = 0.0 PRINT *, " Forming load vector" ! 重力载荷 (本例中V=0,跳过) IF (V .NE. 0.0) THEN CALL PVF(LND,X,Y,P) END IF ! 表面载荷 IF (NPS .GT. 0) THEN CALL PSF(NPS,MPS,ZPS,X,Y,P) END IF ! 集中载荷 (本例中NPJ=0,跳过) IF (NPJ .GT. 0) THEN CALL PCF(NPJ,NPP, PJ,P) END IF ! 保存载荷向量 (可选) ! OPEN(20, FILE='P.TXT') ! DO II = 1, NJ ! WRITE(20,100) 'PX(',II,')=',P(2*II-1),'PY(',II,')=',P(2*II) ! END DO ! 100 FORMAT(1X, 2(A,I3,A,E12.4,3X)) ! CLOSE(20) END SUBROUTINE PFORM ! 重力载荷 (保持不变) SUBROUTINE PVF(LND,X,Y,P) ! 重力 ! ... (与原始代码相同) ... END SUBROUTINE PVF ! 集中力载荷 (保持不变) SUBROUTINE PCF(NPJ,NPP, PJ,P) ! 集中力 ! ... (与原始代码相同) ... END SUBROUTINE PCF ! 表面力载荷 (修复计算) SUBROUTINE PSF(NPS,MPS,ZPS,X,Y,P) ! 分布面力 IMPLICIT NONE INTEGER, INTENT(IN) :: NPS, MPS(:,:) REAL, INTENT(IN) :: ZPS(:,:), X(:), Y(:) REAL, INTENT(INOUT) :: P(:) INTEGER :: NJ, NE, NZ, NPJ, IPS REAL :: E, PR, T, V INTEGER :: I, N1, N2 REAL :: LENGTH, FORCE_X, FORCE_Y COMMON /GLOBAL/ NJ, NE, NZ, NPJ, IPS, E, PR, T, V DO I = 1, NPS IF (I > SIZE(MPS,1)) EXIT N1 = MPS(I,1) N2 = MPS(I,2) ! 添加节点索引检查 IF (N1 < 1 .OR. N1 > NJ .OR. N2 < 1 .OR. N2 > NJ) THEN PRINT *, "WARNING: Invalid node in surface force edge", I CYCLE END IF ! 计算边长 (假设边界是垂直的) LENGTH = ABS(Y(N2) - Y(N1)) ! 计算等效节点力 (均布载荷) FORCE_X = ZPS(I,1) * T * LENGTH / 2.0 ! 每个节点承受一半的力 FORCE_Y = ZPS(I,2) * T * LENGTH / 2.0 ! 施加节点力 (检查索引范围) IF (2*N1-1 <= SIZE(P) .AND. 2*N1 <= SIZE(P)) THEN P(2*N1-1) = P(2*N1-1) + FORCE_X ! X方向 P(2*N1) = P(2*N1) + FORCE_Y ! Y方向 END IF IF (2*N2-1 <= SIZE(P) .AND. 2*N2 <= SIZE(P)) THEN P(2*N2-1) = P(2*N2-1) + FORCE_X ! X方向 P(2*N2) = P(2*N2) + FORCE_Y ! Y方向 END IF END DO END SUBROUTINE PSF ! 约束处理 (修复索引问题) SUBROUTINE RKR(NJ,NZ,JZ,AK,P, DXY) ! 约束处理 IMPLICIT NONE INTEGER, INTENT(IN) :: NJ, NZ, JZ(:,:) REAL, INTENT(IN) :: DXY(:,:) REAL, INTENT(INOUT) :: AK(:,:), P(:) INTEGER :: I, IR, J, II REAL :: LARGE_NUMBER = 1.0E20 PRINT *, " Applying", NZ, "constraints with penalty method" DO I = 1, NZ IF (I > SIZE(JZ,1)) EXIT IR = JZ(I,1) ! 添加节点索引检查 IF (IR < 1 .OR. IR > NJ) THEN PRINT *, "WARNING: Invalid constraint node", IR CYCLE END IF DO J = 2, 3 IF(JZ(I,J) .NE. 0) THEN II = 2*(IR-1) + J-1 ! 添加刚度矩阵索引检查 IF (II < 1 .OR. II > SIZE(AK,1)) THEN PRINT *, "WARNING: Constraint index out of bounds", II CYCLE END IF AK(II,II) = AK(II,II) * LARGE_NUMBER P(II) = AK(II,II) * DXY(I,J-1) END IF END DO END DO END SUBROUTINE RKR ! 高斯消去法求解 (修复索引问题) SUBROUTINE AGAUS(A,B,N,X) ! 高斯消去法 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, INTENT(INOUT) :: A(:,:), B(:) REAL, INTENT(OUT) :: X(:) REAL, ALLOCATABLE :: AM(:,:) REAL :: T, FACTOR INTEGER :: K, I, J ! 检查矩阵维度 IF (SIZE(A,1) < N .OR. SIZE(A,2) < N .OR. SIZE(B) < N .OR. SIZE(X) < N) THEN PRINT *, "ERROR: Matrix dimensions mismatch in AGAUS" STOP END IF ALLOCATE(AM(N,N)) AM = A(1:N,1:N) ! 保存原始矩阵 PRINT *, " Solving system with Gaussian elimination" ! 向前消元 DO K = 1, N-1 ! 检查主元是否为零 IF (ABS(A(K,K)) < 1e-10) THEN PRINT *, "Zero pivot at ", K ! 尝试寻找非零主元 DO I = K+1, N IF (ABS(A(I,K)) > 1e-10) THEN ! 交换 DO J = K, N T = A(K,J) A(K,J) = A(I,J) A(I,J) = T END DO T = B(K) B(K) = B(I) B(I) = T EXIT END IF END DO ! 再次检查 IF (ABS(A(K,K)) < 1e-10) THEN PRINT *, "Zero pivot cannot be resolved at ", K STOP END IF END IF DO I = K+1, N FACTOR = A(I,K)/A(K,K) DO J = K+1, N A(I,J) = A(I,J) - FACTOR*A(K,J) END DO B(I) = B(I) - FACTOR*B(K) END DO END DO ! 回代 X(N) = B(N)/A(N,N) DO I = N-1, 1, -1 T = 0.0 DO J = I+1, N T = T + A(I,J)*X(J) END DO X(I) = (B(I) - T)/A(I,I) END DO A(1:N,1:N) = AM ! 恢复原始矩阵 DEALLOCATE(AM) END SUBROUTINE AGAUS ! 应力计算和输出 (修复索引问题) SUBROUTINE MADE(LND,D,X,Y,U) ! 应力计算 IMPLICIT NONE INTEGER, INTENT(IN) :: LND(:,:) REAL, INTENT(IN) :: D(3,3), X(:), Y(:), U(:) INTEGER :: NJ, NE, NZ, NPJ, IPS REAL :: E, PR, T, V INTEGER :: IE, I, J, IH, IW REAL :: AE, ST(3), UE(6), STX, STY, TXY, AST, RST, STMA, STMI REAL :: B(3,6), S(3,6) COMMON /GLOBAL/ NJ, NE, NZ, NPJ, IPS, E, PR, T, V OPEN(8, FILE='result.txt') WRITE(8,*) 'NODE DISPLACEMENTS AND ELEMENT STRESSES' WRITE(8,10) 10 FORMAT(1X, 'ELEMENT', 6X, 'STX', 14X, 'STY', 14X, 'TXY', 14X, 'ST_MAX', 11X, 'ST_MIN') DO IE = 1, NE IF (IE > SIZE(LND,1)) CYCLE ! 检查单元节点索引 IF (ANY(LND(IE,:) < 1) .OR. ANY(LND(IE,:) > NJ)) THEN WRITE(8,*) 'Element', IE, 'has invalid node indices' CYCLE END IF CALL MA(IE,LND,X,Y,AE) CALL MB(IE,LND,X,Y,AE,B) CALL MS(D, B, S) ! 提取单元位移 UE = 0.0 DO I = 1, 3 DO J = 1, 2 IH = 2*(I-1)+J IW = 2*(LND(IE,I)-1)+J ! 检查位移索引范围 IF (IW >= 1 .AND. IW <= SIZE(U)) THEN UE(IH) = U(IW) ELSE UE(IH) = 0.0 END IF END DO END DO ! 计算应力 ST = 0.0 DO I = 1, 3 DO J = 1, 6 ST(I) = ST(I) + S(I,J)*UE(J) END DO END DO ! 计算主应力 STX = ST(1) STY = ST(2) TXY = ST(3) AST = (STX+STY)*0.5 RST = SQRT(0.25*(STX-STY)**2 + TXY*TXY) STMA = AST + RST STMI = AST - RST WRITE(8,60) IE, STX, STY, TXY, STMA, STMI END DO 60 FORMAT(1X,I5,5E16.6) WRITE(8,70) 70 FORMAT(/,1X,'NODE DISPLACEMENTS',/,' NODE',10X,'UX',14X,'UY') DO I = 1, NJ ! 检查位移索引范围 IF (2*I-1 <= SIZE(U) .AND. 2*I <= SIZE(U)) THEN WRITE(8,80) I, U(2*I-1), U(2*I) ELSE WRITE(8,80) I, 0.0, 0.0 END IF END DO 80 FORMAT(1X,I5,2E16.6) CLOSE(8) PRINT *, " Results saved to result.txt" END SUBROUTINE MADE END PROGRAM PLATE ``` ### 关键修复说明 1. **网格生成问题修复**: - 添加了节点映射表`NODE_MAP`,正确处理孔洞区域 - 使用二维数组跟踪有效节点位置 - 正确处理孔边界的不完整单元 ```fortran INTEGER, ALLOCATABLE :: NODE_MAP(:,:) ! 节点映射表 DO J = 1, NY DO I = 1, NX IF (.NOT. IN_HOLE) THEN NJ = NJ + 1 NODE_MAP(I, J) = NJ END IF END DO END DO ``` 2. **索引越界保护**: - 在所有关键位置添加了数组索引检查 - 防止无效节点导致程序崩溃 ```fortran ! 在AKFORM中 IF (NODE_I < 1 .OR. NODE_I > NJ) THEN PRINT *, "ERROR: Invalid node index in element", IE, "node", I CYCLE END IF ``` 3. **约束处理改进**: - 使用容差比较浮点数,避免精度问题 - 添加约束数量验证 ```fortran NZ = COUNT(ABS(X - 0.0) < 1E-3) ! 使用容差 IF (COUNT_CONSTRAINT < NZ) THEN PRINT *, " WARNING: Only found", COUNT_CONSTRAINT, "constraint nodes" NZ = COUNT_CONSTRAINT ! 更新实际约束数 END IF ``` 4. **面力计算修正**: - 修正了等效节点力的计算方法 - 添加了节点索引检查 ```fortran ! 计算等效节点力 LENGTH = ABS(Y(N2) - Y(N1)) FORCE_X = ZPS(I,1) * T * LENGTH / 2.0 ``` 5. **高斯消去法增强**: - 添加主元为零时的交换处理 - 添加矩阵维度检查 ```fortran ! 检查主元是否为零 IF (ABS(A(K,K)) < 1e-10) THEN ! 尝试寻找非零主元并交换 DO I = K+1, N IF (ABS(A(I,K)) > 1e-10) THEN ! 交换 ... END IF END DO END IF ``` 6. **应力计算保护**: - 添加单元节点有效性检查 - 防止无效位移索引导致错误 ```fortran ! 检查单元节点索引 IF (ANY(LND(IE,:) < 1) .OR. ANY(LND(IE,:) > NJ)) THEN WRITE(8,*) 'Element', IE, 'has invalid node indices' CYCLE END IF ``` ### 验证方法 1. **小规模测试**: ```fortran ! 临时修改网格参数 INTEGER, PARAMETER :: NX = 5, NY = 3 ! 小网格 ``` 2. **结果验证**: - 检查孔边应力集中系数 ≈ 3.0 - 验证右端位移与理论值一致 - 检查对称性(Y=0.5m处应力分布对称) 3. **调试输出**: ```fortran ! 在关键子程序添加调试输出 PRINT *, "Processing element", IE, "with nodes:", LND(IE,:) ``` ### 预期结果 1. **变形特征**: - 左端固定,位移接近零 - 右端均匀拉伸,X向位移最大 - 孔洞区域变形复杂,呈现椭圆形 2. **应力分布**: - 孔边水平直径两端应力最大 - 应力集中系数 ≈ 3.0 - 远离孔洞区域应力趋于均匀 3. **数值验证**: - 最大位移:约 2.5e-5 m - 最大应力:约 60,000 Pa (3 × 20,000 Pa) ### 使用建议 1. **网格细化**: - 对于精确的应力分析,建议使用更细的网格 - 孔边网格尺寸 ≤ 0.05m 2. **参数调整**: ```fortran ! 在GENERATE_MESH_WITH_HOLE中修改 INTEGER, PARAMETER :: NX = 80, NY = 30 ! 更高分辨率 ``` 3. **结果可视化**: ```python # Python结果分析脚本 import numpy as np import matplotlib.pyplot as plt # 读取结果文件 data = np.loadtxt('result.txt', skiprows=10) plt.figure(figsize=(12, 6)) plt.subplot(121) plt.scatter(data[:, 0], data[:, 1], c=data[:, 4], cmap='jet') plt.colorbar(label='Max Principal Stress (Pa)') plt.title('Stress Distribution') plt.subplot(122) plt.scatter(data[:, 0], data[:, 1], c=data[:, 5], cmap='viridis') plt.colorbar(label='Displacement Magnitude (m)') plt.title('Displacement Field') plt.tight_layout() plt.savefig('plate_results.png') ``` ---
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值