     !     This is part of the netCDF package.
     !     Copyright 2006 University Corporation for Atmospheric
     !     Research/Unidata.
     !     See COPYRIGHT file for conditions of use.

     !     This is a very simple example which writes a 2D array of
     !     sample
     !     data. To handle this in netCDF we create two shared
     !     dimensions,
     !     "x" and "y", and a netCDF variable, called "data". It uses
     !     parallel I/O to write the file from all processors at the
     !     same
     !     time.

     !     This example demonstrates the netCDF Fortran 90 API. This is
     !     part
     !     of the netCDF tutorial, which can be found at:
     !     http://www.unidata.ucar.edu/software/netcdf/docs/netcdf-tutorial

     !     Full documentation of the netCDF Fortran 90 API can be found
     !     at:
     !     http://www.unidata.ucar.edu/software/netcdf/docs/netcdf-f90

     !     $Id: simple_xy_par_wr.f90,v 1.2 2009/03/12 18:29:48 ed Exp $

     program simple_xy_par_wr
       use netcdf
       implicit none
       include 'mpif.h'

       ! This is the name of the data file we will create.
       character (len = *), parameter :: FILE_NAME = "simple_xy_par.nc"

       ! We are writing 2D data.
       integer, parameter :: NDIMS = 2

       ! When we create netCDF files, variables and dimensions, we get
       ! back
       ! an ID for each one.
       integer :: ncid, varid, dimids(NDIMS)
       integer :: x_dimid, y_dimid

       ! These will tell where in the data file this processor should
       ! write.
       integer :: start(NDIMS), count(NDIMS)

       ! This is the data array we will write. It will just be filled
       ! with
       ! the rank of this processor.
       integer, allocatable :: data_out(:)

       ! MPI stuff: number of processors, rank of this processor, and
       ! error
       ! code.
       integer :: p, my_rank, ierr

       ! Loop indexes, and error handling.
       integer :: x, y, stat

       ! Initialize MPI, learn local rank and total number of
       ! processors.
       call MPI_Init(ierr)
       call MPI_Comm_rank(MPI_COMM_WORLD, my_rank, ierr)
       call MPI_Comm_size(MPI_COMM_WORLD, p, ierr)

       ! Create some pretend data. We just need one row.
       allocate(data_out(p), stat = stat)
       if (stat .ne. 0) stop 3
       do x = 1, p
          data_out(x) = my_rank
       end do

       ! Create the netCDF file. The NF90_NETCDF4 flag causes a
       ! HDF5/netCDF-4 file to be created. The comm and info parameters
       ! cause parallel I/O to be enabled.
       call check( nf90_create(FILE_NAME, NF90_NETCDF4, ncid,  &
                   comm=MPI_COMM_WORLD, info=MPI_INFO_NULL) )

       ! Define the dimensions. NetCDF will hand back an ID for
       ! each. Metadata operations must take place on all processors.
       call check( nf90_def_dim(ncid, "x", p, x_dimid) )
       call check( nf90_def_dim(ncid, "y", p, y_dimid) )

       ! The dimids array is used to pass the IDs of the dimensions of
       ! the variables. Note that in fortran arrays are stored in
       ! column-major format.
       dimids =  (/ y_dimid, x_dimid /)

       ! Define the variable. The type of the variable in this case is
       ! NF90_INT (4-byte integer).
       call check( nf90_def_var(ncid, "data", NF90_INT, dimids, varid) )

       ! End define mode. This tells netCDF we are done defining
       ! metadata. This operation is collective and all processors will
       ! write their metadata to disk.
       call check( nf90_enddef(ncid) )

       ! Write the pretend data to the file. Each processor writes one
       ! row.
       start = (/ 1, my_rank + 1/)
       count = (/ p, 1 /)
       call check( nf90_put_var(ncid, varid, data_out, start = start, &
            count = count) )

       ! Close the file. This frees up any internal netCDF resources
       ! associated with the file, and flushes any buffers.
       call check( nf90_close(ncid) )

       ! Free my local memory.
       deallocate(data_out)

       ! MPI library must be shut down.
       call MPI_Finalize(ierr)

       if (my_rank .eq. 0)  then
         print *, "*** SUCCESS writing example file ", FILE_NAME, "! "
       endif

     contains
       subroutine check(status)
         integer, intent ( in) :: status

         if(status /= nf90_noerr) then
           print *, trim(nf90_strerror(status))
           stop 2
         end if
       end subroutine check
     end program simple_xy_par_wr


text:
program rstest
  use flips
  use rand

  ! This program tests flips_resize in different situations
  ! Single precision real version

  implicit none

  type(flips_d) :: s1,s2,s3,s4

  real(dp), dimension(:), allocatable :: amat,amat2,rr,meas2
  real(dp), dimension(:), allocatable :: meas,erro
  logical, dimension(:), allocatable :: re
  real(dp) :: diff, maxdiff
  real(dp) :: ts,te
  logical, dimension(:), allocatable :: remo
  integer :: i

  call set_rand_seed()

  write(*,*) 'Test 1: Expand 1000 unknowns problem by 1000 new unknowns'
  ! ********************************************************************

  write(*,*) '  Memory storage'

  ! Initialize s1
  call flips_init(s1,1000,1,buffersize=50)
  call flips_init(s3,500,1,buffersize=50)

  ! Create data
  allocate(amat(1000*1000),meas(1000),erro(1000),amat2(1500*500),rr(500*500),meas2(500),remo(1000))
  call random_number(amat)
  call random_number(meas)
  !call random_number(erro)
  erro = 1.0

  ! Feed data in s1
  call flips_add(s1,1000,amat,meas,erro)

  call cpu_time(ts)

  remo = .FALSE.

  call flips_resize(s2,s1,newsize=1500,remove=remo)

 ! write(*,*) s2%ymat
  write(*,*) 'ncols',s2%ncols
  write(*,*) 'nrhs',s2%nrhs
  write(*,*) 'nbuf',s2%nbuf
  write(*,*) 'nrows',s2%nrows
  write(*,*) 'bw',s2%bw
  write(*,*) 'bbw',s2%bbw
  write(*,*) 'common',s2%common
  write(*,*) 'rl',s2%rl
  write(*,*) 'zeroth',s2%zeroth
  write(*,*) 'nrotbuf',s2%nrotbuf

  !s2%nrows = 1000

  call cpu_time(te)

  ! Add some new data

  amat2 = 0.0
  call random_number(rr)
  do i = 1,500
     rr(yind(i,i,500)) = 1.0
  end do

  call random_number(meas2)

  do i = 1,500

     amat2(yind(i,1001,1500):yind(i,1500,1500)) = rr(yind(i,1,500):yind(i,500,500))

     !write(*,*) amat2(yind(i,996,1500):yind(i,1005,1500))

  end do

     !write(*,*) amat2(996:1005)

     !call random_number(meas2)

     call flips_add(s2,500,amat2,meas2,erro(1:500))
     call flips_add(s3,500,rr,meas2,erro(1:500))

  call flips_rotate(s1)
  call flips_rotate(s2)
  call flips_rotate(s3)

 ! call flips_resize(s4,s2,newsize=500)

  !write(*,*) s2%ymat(1:1000)-s1%ymat

  ! Solve all problems
  call flips_solve(s1)

  call flips_solve(s2)

  call flips_solve(s3)

 ! call flips_solve(s4)

write(*,*) s1%solmat(1:5)
write(*,*) s2%solmat(1:5)
write(*,*) s2%solmat(1001:1005)
write(*,*) s3%solmat(1:5)
!write(*,*) s4%solmat(1:5)

  ! Compare results
  diff = sum(s2%solmat(1:1000) - s1%solmat)/1000
  maxdiff = maxval(s2%solmat(1:1000) - s1%solmat)
  write(*,*) '    Mean diff:',diff
  write(*,*) '    Max diff:',maxdiff
  write(*,*) '    Resize time:',te-ts,'\n'

  call flips_kill(s1,.FALSE.)
  call flips_kill(s2,.FALSE.)
  call flips_kill(s3,.FALSE.)
 ! call flips_kill(s4,.FALSE.)
  deallocate(amat,amat2,meas,erro,rr,meas2)

end program rstest

max_stars_repo_name: morispaa/flips

max_issues_repo_name: morispaa/flips

max_forks_repo_name: morispaa/flips


text:
*     SB09MD EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2020 NICONET e.V.
*
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        ( NIN = 5, NOUT = 6 )
      INTEGER          NMAX, NCMAX, NBMAX
      PARAMETER        ( NMAX = 20, NCMAX = 20, NBMAX = 20 )
      INTEGER          LDH1, LDH2, LDSS, LDSE, LDPRE
      PARAMETER        ( LDH1 = NCMAX, LDH2 = NCMAX, LDSS = NCMAX,
     $                   LDSE = NCMAX, LDPRE = NCMAX )
*     .. Local Scalars ..
      DOUBLE PRECISION TOL
      INTEGER          I, INFO, J, N, NB, NC
*     .. Local Arrays ..
      DOUBLE PRECISION H1(LDH1,NMAX*NBMAX), H2(LDH2,NMAX*NBMAX),
     $                 PRE(LDPRE,NBMAX), SE(LDSE,NBMAX), SS(LDSS,NBMAX)
*     .. External Subroutines ..
      EXTERNAL         SB09MD
*     .. Executable Statements ..
*
      WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read the data.
      READ ( NIN, FMT = '()' )
      READ ( NIN, FMT = * ) N, NC, NB, TOL
      IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
         WRITE ( NOUT, FMT = 99993 ) N
      ELSE IF ( NB.LT.0 .OR. NB.GT.NBMAX ) THEN
         WRITE ( NOUT, FMT = 99992 ) NB
      ELSE IF ( NC.LT.0 .OR. NC.GT.NCMAX ) THEN
         WRITE ( NOUT, FMT = 99991 ) NC
      ELSE
         READ ( NIN, FMT = * ) ( ( H1(I,J), I = 1,NC ), J = 1,N*NB )
         READ ( NIN, FMT = * ) ( ( H2(I,J), I = 1,NC ), J = 1,N*NB )
*        Compare the given sequences and evaluate their closeness.
         CALL SB09MD( N, NC, NB, H1, LDH1, H2, LDH2, SS, LDSS, SE, LDSE,
     $                PRE, LDPRE, TOL, INFO )
*
         IF ( INFO.NE.0 ) THEN
            WRITE ( NOUT, FMT = 99998 ) INFO
         ELSE
            WRITE ( NOUT, FMT = 99997 )
            DO 20 I = 1, NC
               WRITE ( NOUT, FMT = 99996 ) ( SS(I,J), J = 1,NB )
   20       CONTINUE
            WRITE ( NOUT, FMT = 99995 )
            DO 40 I = 1, NC
               WRITE ( NOUT, FMT = 99996 ) ( SE(I,J), J = 1,NB )
   40       CONTINUE
            WRITE ( NOUT, FMT = 99994 )
            DO 60 I = 1, NC
               WRITE ( NOUT, FMT = 99996 ) ( PRE(I,J), J = 1,NB )
   60       CONTINUE
         END IF
      END IF
      STOP
*
99999 FORMAT (' SB09MD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from SB09MD = ',I2)
99997 FORMAT (' The sum-of-squares matrix SS is ')
99996 FORMAT (20(1X,F8.4))
99995 FORMAT (/' The quadratic error matrix SE is ')
99994 FORMAT (/' The percentage relative error matrix PRE is ')
99993 FORMAT (/' N is out of range.',/' N = ',I5)
99992 FORMAT (/' NB is out of range.',/' NB = ',I5)
99991 FORMAT (/' NC is out of range.',/' NC = ',I5)
      END

ext: f

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: bnavigator/SLICOT-Reference

max_forks_repo_name: bnavigator/SLICOT-Reference


text:
    module atp_util_mod
    use sim_tensor_mod
    use sim_util_mod
    use atp_element_mod
    use output_mod
    implicit none

    private

    public  ::  solve_incr      ! Attempt to solve one increment
    public  ::  apply_load      ! Apply the external loads
    public  ::  apply_bc        ! Apply the boundary conditions (prescribed displacements)
    public  ::  get_result      ! Get the results from the particular increment
    public  ::  gen_free_dofs   ! Obtain the free degrees of freedom for particular step
    public  ::  setup_solve_incr    ! Initiate variables for solve_incr
    public  ::  get_equilibrium_didnt_converge_count ! As name implies

    ! Parameters
    double precision, parameter :: pi = acos(-1.d0)  !Define pi

    ! Simulation variables
    integer         :: equilibrium_didnt_converge_count     ! Count number of times equilibrium did not converge
    logical         :: full_output                          ! Should the info be output every time no convergence?

    contains
    !==========================================================================================================
    !Iteration subroutines

subroutine setup_solve_incr(set_full_output)
implicit none
    logical set_full_output
    full_output = set_full_output
    equilibrium_didnt_converge_count = 0
end subroutine

!Solve one increment
subroutine solve_incr(Rpos, H0, load, disp, temp, dtemp, time, dt, free_dofs, known_dofs, gp_F, gp_s, gp_sv, gp_strain, u, du, &
    iter, niter, lconv, pnewdt, kinc, kstep, props, cmname, umat_address, nlgeom, iter_err_norm, denergy)
    use types_mod
    use umat_mod
    implicit none
        !Input/Output variables
        double precision, intent(in)    :: Rpos(:,:), H0   !Element nodal coordinates, dim = (nnod, nel) (nnod = number of nodes in element)
        double precision, intent(inout) :: load(:), disp(:) !External loads and displacements
        double precision, intent(in)    :: temp, dtemp, time(2), dt ! Time and temperature
        double precision, intent(inout) :: gp_F(:,:,:), gp_s(:,:,:), gp_sv(:,:,:), gp_strain(:,:,:)  !Values at gauss points (dim(*),ngp,nel)
        double precision, intent(inout) :: u(:), du(:)  !Displacements (uz, phi, ur1,ur2,ur3,...,urN)
        type(iter_typ)                  :: iter         ! Iteration settings
        double precision, intent(out)   :: pnewdt       ! Variable used to request shorter time step
        double precision, intent(in)    :: props(:)     ! Material parameters
        integer, intent(out)            :: niter        ! Number of iterations used
        integer, intent(in)             :: kinc, kstep  ! Current increment and step
        integer, intent(in)             :: free_dofs(:) ! List over the non-prescribed (i.e. free) degrees of freedom
        integer, intent(in)             :: known_dofs(:)! List over the degrees of freedom with known load
        logical, intent(out)            :: lconv        ! Local convergence bool
        character(len=80), intent(in)   :: cmname       ! Material name sent to umat
        procedure(umat_template),pointer:: umat_address ! Addresss to umat subroutine
        logical, intent(in)             :: nlgeom       ! Bool to determine if nonlinear geometry effects should be accounted for
        double precision, intent(in)    :: iter_err_norm(:)	!Scaling to obtain averaged stresses from the residuals
        double precision, intent(out)   :: denergy(4)   ! sse, spd, scd, rpl increment, integrated over the body (i.e. total energy, not specific)

        !Internal variables
        double precision :: K(size(u), size(u)), R(size(u))
        double precision :: error, error_old, alpha
        double precision :: Ke(size(Rpos,1)+2, size(Rpos,1)+2), Re(size(Rpos,1)+2)
        integer          :: nel, nnod, ndof_tot, k1, k2, dof(size(Rpos,1)+2)
        double precision :: gp_F0(9, size(gp_F, dim=2), size(gp_F, dim=3))
        double precision :: gp_s0(6, size(gp_s, dim=2), size(gp_s, dim=3))
        double precision :: gp_strain0(6, size(gp_strain,dim=2), size(gp_strain,dim=3))
        double precision :: u_old(size(u)), ue(size(Rpos,1)+2)
        double precision :: Rscale      ! Scale factor for the equilibrium equations
        double precision, allocatable :: gp_sv0(:,:,:)
        double precision :: elem_denergy(4)         ! sse, spd, scd, rpl increment, integrated over element

        allocate(gp_sv0(size(gp_sv,dim=1), size(gp_sv,dim=2), size(gp_sv,dim=3)))

        nnod= size(Rpos,1) !Number of nodes per element
        nel = size(Rpos,2) !Number of elements
        ndof_tot = 2 + 1+(nnod-1)*nel

        lconv = .true.

        dof(1:2) = (/1,2/) ! All elements have the axial and torsion dofs

        error_old = huge(1.d0)  ! Ensure that error is decreasing the first time, to avoid line search.
        alpha = 1.d0            ! Initial alpha line search factor

        u_old = u   !Save old u to make it possible to calculate total du as input to the next step's guess
        u = u + du  !du contains correct displacements for prescribed displacement increment

        gp_F0 = gp_F
        gp_s0 = gp_s
        gp_sv0= gp_sv
        gp_strain0 = gp_strain
        pnewdt = 1.d0
        niter = 0
        do while(.true.)
            ! Step 0: Reset values for each iteration attempt
            gp_F  = gp_F0
            gp_s  = gp_s0
            gp_sv = gp_sv0
            gp_strain = gp_strain0
            denergy = 0.d0

            R       = 0.d0
            Rscale  = 0.d0
            K       = 0.d0

            ! Step 1: Assemble system
            do k1=1,nel
                do k2=1,nnod
                    dof(k2+2) = (k2+2) + (nnod-1)*(k1-1)
                enddo

                ! Get element displacements
                ue = u(dof)

                ! Call element routine
                if (nlgeom) then
                    call element_nlgeom( Ke, Re, ue, gp_s(:,:,k1), gp_sv(:,:,k1), gp_F(:,:,k1), Rpos(:,k1), H0, &
                    kinc, kstep, k1, pnewdt, props, cmname, dtemp, temp, time, dt, umat_address, elem_denergy)
                else
                    call element_lingeom( Ke, Re, ue, gp_strain(:,:,k1), gp_s(:,:,k1), gp_sv(:,:,k1), gp_F(:,:,k1), Rpos(:,k1), H0, &
                    kinc, kstep, k1, pnewdt, props, cmname, dtemp, temp, time, dt, umat_address, elem_denergy)
                endif

                if (pnewdt<1.d0) then
                    lconv = .false.
                    exit
                endif

                K(dof,dof)  = K(dof, dof) + Ke
                R(dof)      = R(dof)      + Re
                Rscale      = Rscale + sqrt(sum((Re*iter_err_norm(dof))**2))	! Add the average stress contributions to Rscale
                denergy = denergy + elem_denergy
            enddo
            if (.not.lconv) then
                exit
            endif

            ! Apply external loads (include potential effect on stiffness if nlgeom=true)
            call apply_load(R, K, load, H0, Rpos, u, known_dofs, nlgeom)

            ! Step 2: Check residual
            if (free_dofs(1)/=-1) then
                error = sqrt(sum((R(known_dofs)*iter_err_norm(known_dofs))**2))/max(Rscale,1.d0)
            !   error = sqrt(sum(R(known_dofs)**2))
            else
                error = 0.d0
            endif

            ! Step 2: Solve increment update, depending on error
            if (error<iter%tol) then
                lconv = .true.
                exit
            elseif ((error>error_old).and.(niter>iter%ls_start)) then
                u = u - du
                alpha = alpha*iter%ls_alpha_factor
                du = alpha*du
            else !tol<error<error_old
                alpha = 1.d0
                call solve_mateq(K, R, du, lconv, free_dofs, known_dofs)

                if (.not.lconv) then
                    pnewdt = 0.5d0
                    call write_output('Could not solve matrix equation, attempting a smaller time increment', 'status', 'sim:atp')
                    lconv = .false.
                    exit
                endif
            endif

            ! Step 3: Update solution
            u = u + du
            niter = niter + 1
            error_old = error

            ! Step 4: Check for max number of iterations
            if (niter>iter%max) then
                lconv = .false.
                pnewdt = 0.75d0
                if (full_output) then
                    call write_output('Could not find equilibrum in '//int2str(iter%max)//' attempts, using smaller time step', 'status', 'sim:atp')
                else
                    equilibrium_didnt_converge_count = equilibrium_didnt_converge_count + 1
                endif
                exit
            endif
        enddo

        du = u - u_old  !Output total displacement increment

        ! Insert the previously unknown values into load/disp vectors
        call get_result(R, load, disp, H0, Rpos, u, free_dofs, known_dofs, nlgeom)

end subroutine

!Load and boundary conditions
subroutine apply_load(R, K, load, H0, Rpos, u, known_dofs, nlgeom)
    implicit none
    double precision    :: R(:), K(:,:), load(:), H0, Rpos(:,:), u(:)
    integer             :: known_dofs(:), n
    logical             :: nlgeom
    double precision    :: h, ri, ro

    n    = size(u)

    h = H0                                  ! Height
    ri = Rpos(1,1)                          ! Inner radius
    ro = Rpos(size(Rpos,1),size(Rpos,2))    ! Outer radius

    if (nlgeom) then    ! get current:
        h  = H0 + u(1)  ! height,
        ri = ri + u(3)  ! inner radius and
        ro = ro + u(n)  ! outer radius
    endif

    ! Axial force
    if (any(known_dofs==1)) then
        R(1) = R(1) - load(1)
    endif

    ! Torque
    if (any(known_dofs==2)) then
        R(2) = R(2) - load(2)
    endif

    ! Internal pressure
    if (any(known_dofs==3)) then
        R(3) = R(3) - load(3)*(2.d0*ri*pi)*h
        if (nlgeom) then
            K(3,1)  = K(3,1) - load(3)*(2.d0*ri*pi)
            K(3,3)  = K(3,3) - load(3)*(2.d0*pi)*h
        endif
    endif

    ! External pressure
    if (any(known_dofs==n)) then
        R(n) = R(n) + load(4)*(2.d0*ro*pi)*h
        if (nlgeom) then
            K(n,1)  = K(n,1) + load(4)*(2.d0*ro*pi)
            K(n,n)  = K(n,n) + load(4)*(2.d0*pi)*h
        endif
    endif

end subroutine

subroutine apply_bc(disp, disp_new, du, ctrl, disp_conv)
    implicit none
    double precision, intent(in)    :: disp(:), disp_new(:), disp_conv(:)
    integer, intent(in)             :: ctrl(:)
    double precision                :: du(:)
    integer                         :: ctrl_dofs(size(ctrl)), k1

    ctrl_dofs = (/1, 2, 3, size(du)/)

    do k1=1,4
        if (ctrl(k1)>0) then !Displacement control
            du(ctrl_dofs(k1)) = (disp_new(k1)-disp(k1))*disp_conv(k1)
        endif
    enddo

end subroutine

subroutine get_result(R, load, disp, H0, Rpos, u, free_dofs, known_dofs, nlgeom)
    implicit none
    double precision    :: R(:), load(:), disp(:), H0, Rpos(:,:), u(:)
    integer             :: free_dofs(:), known_dofs(:)
    logical             :: nlgeom
    integer             :: n, pos(4), k1
    double precision    :: h, ri, ro, disp_conv(4)

    n    = size(u)
    pos = (/1,2,3,n/)
    h = H0                                  ! Height
    ri = Rpos(1,1)                          ! Inner radius
    ro = Rpos(size(Rpos,1),size(Rpos,2))    ! Outer radius

    ! Note that this is not the same disp_conv as the globally defined (if ri<1.d-12)
    disp_conv = (/H0, 1.d0, ri, ro/)        ! Must be before nlgeom changes below.
    if (ri<1.d-12) then  ! If inner radius = 0, then strain is infinity. Hence the output is the inner displacement which should be zero!
        disp_conv(3) = huge(1.d0)
    endif

    if (nlgeom) then    ! get current:
        h  = H0 + u(1)  ! height,
        ri = ri + u(3)  ! inner radius and
        ro = ro + u(n)  ! outer radius
    endif

    ! Update outputs
    do k1 = 1,4
        if (any(free_dofs==pos(k1))) then       ! Displacement to be determined
            disp(k1) = u(pos(k1))/disp_conv(k1) ! Obtain resulting strains
        endif
        if (.not.any(known_dofs==pos(k1))) then ! Force to be determined
            if (pos(k1)<=2) then                ! Axial force or torque
                load(k1) = R(pos(k1))           ! Obtain resulting load for axial force or torque
            elseif (pos(k1)==3) then            ! Inner pressure
                if (ri<1.d-12) then             ! If inner radius is equal to zero (i.e. very small)
                    load(3) = 0.d0              ! Set load to zero
                else
                    load(3) = R(pos(k1))/(2*pi*ri*h)
                endif
            else                                ! Outer pressure
                load(4) = -R(pos(k1))/(2*pi*ro*h)
            endif
        endif
    enddo

end subroutine

!Generate free_dofs
subroutine gen_free_dofs(free_dofs, known_dofs, ndof_tot, ctrl, intpres_extstrn)
    implicit none
        integer, allocatable :: free_dofs(:)            ! The free degrees of freedom
        integer, allocatable :: known_dofs(:)           ! The degrees of freedom with known load
        integer              :: ctrl(4)                 ! The control parameters (axial, torsion, inside, outside)
        logical              :: intpres_extstrn         ! If internal pressure and external strain used with strain control
        integer              :: ctrl_dofs(4)            ! Degree of freedom for dofs directly controlled by ctrl
        integer              :: ndof_tot, nfree
        integer              :: k1, k2, k3

        nfree = ndof_tot-count(ctrl>0)  ! Get number of free dofs

        if (nfree==0) then  ! If no free degrees of freedom, we handle this special case separately!
            allocate(free_dofs(1))
            free_dofs = -1
            allocate(known_dofs, source=free_dofs)
        else
            allocate(free_dofs(nfree))
            ctrl_dofs = (/1, 2, 3, ndof_tot/)
            ! Find free dofs amongst the first specified in ctrl
            k2 = 1
            do k1 = 1,3
                if (ctrl(k1)<0) then ! Force controlled
                    free_dofs(k2) = ctrl_dofs(k1)
                    k2 = k2 + 1
                endif
            enddo

            if (k2<=nfree) then  ! There are additional dofs
                k3 = 4  ! Three dofs (axial, tors, inner are before the additional dofs)
                do k1=k2,nfree
                    free_dofs(k1) = k3
                    k3 = k3 + 1
                enddo
            endif

            ! We want sorted dofs, hence if ctrl(4)<0 then this dof should be put at the end
            ! Then the last is overwritten, but it is actually ndof_tot so the code below is not needed:
            !if (ctrl(4)<0) then
            !    free_dofs(nfree) = ndof_tot
            !endif
            ! Normally, the known_dofs (load known) corresponds to free_dofs (disp unknown)
            allocate(known_dofs, source=free_dofs)
            ! In the special case of known external pressure and external strain,
            ! with unknown internal disp and internal pressure, they are different:
            if (intpres_extstrn) then
                k1 = count(ctrl(1:2)<0)+1
                if (free_dofs(nfree)==ndof_tot) then
                    ! For special intpres_extstrn mode the external strain should be controlled
                    call write_output('For the special intpres_extstrn mode the external strain must be controlled', 'error', 'sim:atp')
                elseif (free_dofs(k1)/=3) then
                    call write_output('For the special intpres_extstrn mode the internal pressure should be set as controlled', 'error', 'sim:atp')
                endif

                known_dofs(k1:) = free_dofs(k1:)+1
            endif

        endif

end subroutine

function get_equilibrium_didnt_converge_count() result (equilibrium_didnt_converge_count_out)
implicit none
    integer equilibrium_didnt_converge_count_out
    equilibrium_didnt_converge_count_out = equilibrium_didnt_converge_count
end function

end module

max_stars_repo_name: KnutAM/matmodfit

max_issues_repo_name: KnutAM/matmodfit

max_forks_repo_name: KnutAM/matmodfit