Fortran 实现整体自由度序号指示矩阵：二维平面节点约束示例

下面是整体自由度序号指示矩阵的 Fortran 语言实现：

integer :: ndof, nnode, ncons, i, j, k, ii, jj, kk
integer, dimension(:), allocatable :: node_dof, node_cons
integer, dimension(:,:), allocatable :: index_matrix

! 定义节点自由度和约束情况
ndof = 2 ! 每个节点有两个自由度
nnode = 6 ! 共有6个节点
ncons = 4 ! 共有4个约束点
allocate(node_dof(ndof*nnode))
allocate(node_cons(ncons))
node_dof = reshape([(i-1)*ndof+1:i*ndof, i=1,nnode])
node_cons = [4,5,6,3] ! 约束点的编号

! 初始化指示矩阵
allocate(index_matrix(ndof*nnode-2*ncons, nnode*ndof))
index_matrix = 0

! 构造指示矩阵
k = 0
do i = 1, nnode
  ! 判断节点是否约束
  if (any(node_cons == i)) then
    continue
  endif
  k = k + 1
  do j = 1, ndof
    index_matrix(k, node_dof(j,i)) = 1
  enddo
enddo

! 输出指示矩阵
do i = 1, size(index_matrix,1)
  do j = 1, size(index_matrix,2)
    write(*,*) index_matrix(i,j),
  enddo
  write(*,*)  
enddo

其中，ndof表示每个节点的自由度个数，本例中为2；nnode表示节点个数，本例中为6；ncons表示约束点个数，本例中为4；node_dof是一个长度为ndof*nnode的一维数组，表示每个节点的自由度编号；node_cons是一个长度为ncons的一维数组，表示约束点的编号；index_matrix是一个大小为(ndof*nnode-2*ncons) x (nnode*ndof)的二维数组，表示指示矩阵。

算法的实现过程如下：

定义节点自由度和约束情况。
初始化指示矩阵。
对于每个未约束的节点，遍历该节点的自由度，将指示矩阵中对应位置的值设为1。
输出指示矩阵。

需要注意的是，指示矩阵的大小为(ndof*nnode-2*ncons) x (nnode*ndof)，其中-2*ncons表示约束点的自由度不需要纳入指示矩阵中。