SetupInterfaces_MultiSetModel3D_t Subroutine

  subroutine SetupInterfaces_MultiSetModel3D_t(this)
    !! Finish setup once every set model has been allocated and initialized.
    !!
    !! Checks that every set model exists, was initialized on its own set's mesh and
    !! geometry, and has gradients disabled. Builds the interface map, then registers a no-op hyperbolic
    !! boundary condition for SELF_BC_SET_INTERFACE on every set model and maps its
    !! boundary conditions again, so interface faces are not counted as unmapped
    !! boundaries.
    !!
    !! Collective: MapBoundaryConditions reduces over each set's communicator, and the
    !! interface map and face count use the global one. Every rank must call it.
    implicit none
    class(MultiSetModel3D_t),intent(inout) :: this
    ! Local
    procedure(SELF_bcMethod),pointer :: bcfunc
    integer :: k,iEl,s,nLocal,ierror
    integer,allocatable :: nvarOfSet(:)

    do k = 1,this%nSets
      if(.not. allocated(this%sets(k)%model)) then
        print*,__FILE__," : set ",k," has no model allocated"
        error stop "MultiSetModel3D : every set needs a model before SetupInterfaces"
      endif
      if(.not. associated(this%sets(k)%model%mesh,this%sets(k)%mesh)) then
        print*,__FILE__," : the model of set ",k," was not initialized on that set's mesh"
        error stop "MultiSetModel3D : initialize sets(k)%model on sets(k)%mesh"
      endif
      if(.not. associated(this%sets(k)%model%geometry,this%sets(k)%geometry)) then
        print*,__FILE__," : the model of set ",k," was not initialized on that set's geometry"
        error stop "MultiSetModel3D : initialize sets(k)%model on sets(k)%geometry"
      endif
      if(this%sets(k)%model%gradient_enabled) then
        print*,__FILE__," : the model of set ",k," has gradient_enabled"
        error stop "MultiSetModel3D : solution gradients are not supported"
      endif
    enddo

    nLocal = 0
    do k = 1,this%nSets
      do iEl = 1,this%sets(k)%mesh%nElem
        do s = 1,6
          if(this%sets(k)%mesh%sideInfo(3,s,iEl) == 0 .and. &
             this%sets(k)%mesh%sideInfo(5,s,iEl) == SELF_BC_SET_INTERFACE) then
            nLocal = nLocal+1
          endif
        enddo
      enddo
    enddo
    call mpi_allreduce(nLocal,this%nInterfaceFaces,1,MPI_INTEGER,MPI_SUM, &
                       this%globalMesh%decomp%mpiComm,ierror)

    allocate(nvarOfSet(1:this%nSets))
    do k = 1,this%nSets
      nvarOfSet(k) = this%sets(k)%model%nvar
    enddo
    call this%iface%BuildMap(this%globalMesh,this%setOfElem,this%setGlobalId, &
                             this%setOffsetThisRank,nvarOfSet, &
                             this%sets(1)%model%solution%interp%N)
    deallocate(nvarOfSet)
    ! The sub-meshes tag a face as an interface exactly when the global mesh puts the
    ! two elements in different sets, so both counts must agree.
    if(this%iface%nFaces /= nLocal) then
      error stop "MultiSetModel3D : interface map disagrees with the sub-mesh tags"
    endif

    allocate(this%traces(1:this%nSets))
    do k = 1,this%nSets
      this%traces(k)%trace => this%sets(k)%model%solution%boundary
    enddo

    ! Size each sub-mesh's request scratch for its own halo before stepping: a send and a
    ! receive per rank-remote side (or mortar sub-side) per variable, as the exchanges
    ! reserve. The sub-meshes are built with a minimal scratch.
    do k = 1,this%nSets
      call ReserveSetMessages(this%sets(k)%mesh,this%sets(k)%model%nvar)
    enddo

    bcfunc => hbc3d_SetInterface_MultiSetModel3D
    do k = 1,this%nSets
      call this%sets(k)%model%hyperbolicBCs%RegisterBoundaryCondition( &
        SELF_BC_SET_INTERFACE,"set_interface",bcfunc)
      call this%sets(k)%model%MapBoundaryConditions()
    enddo

    this%interfacesReady = .true.

  endsubroutine SetupInterfaces_MultiSetModel3D_t