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
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值