sBase_init Subroutine

private subroutine sBase_init(sf, deg_in, continuity_in, grid_in, degGP_in)

Uses

  • proc~~sbase_init~~UsesGraph proc~sbase_init t_sBase%sBase_init module~modgvec_basis1d MODgvec_Basis1D proc~sbase_init->module~modgvec_basis1d module~modgvec_globals MODgvec_Globals proc~sbase_init->module~modgvec_globals module~modgvec_linalg MODgvec_LinAlg proc~sbase_init->module~modgvec_linalg module~sll_m_boundary_condition_descriptors sll_m_boundary_condition_descriptors proc~sbase_init->module~sll_m_boundary_condition_descriptors module~sll_m_bsplines sll_m_bsplines proc~sbase_init->module~sll_m_bsplines module~sll_m_spline_interpolator_1d sll_m_spline_interpolator_1d proc~sbase_init->module~sll_m_spline_interpolator_1d module~sll_m_spline_matrix sll_m_spline_matrix proc~sbase_init->module~sll_m_spline_matrix module~modgvec_basis1d->module~modgvec_globals iso_fortran_env iso_fortran_env module~modgvec_globals->iso_fortran_env module~modgvec_linalg->module~modgvec_globals module~sll_m_bsplines_base sll_m_bsplines_base module~sll_m_bsplines->module~sll_m_bsplines_base module~sll_m_bsplines_non_uniform sll_m_bsplines_non_uniform module~sll_m_bsplines->module~sll_m_bsplines_non_uniform module~sll_m_bsplines_uniform sll_m_bsplines_uniform module~sll_m_bsplines->module~sll_m_bsplines_uniform module~sll_m_working_precision sll_m_working_precision module~sll_m_bsplines->module~sll_m_working_precision module~sll_m_spline_interpolator_1d->module~sll_m_boundary_condition_descriptors module~sll_m_spline_interpolator_1d->module~sll_m_spline_matrix module~sll_m_spline_interpolator_1d->module~sll_m_bsplines_base module~sll_m_spline_1d sll_m_spline_1d module~sll_m_spline_interpolator_1d->module~sll_m_spline_1d module~sll_m_spline_interpolator_1d->module~sll_m_working_precision module~sll_m_spline_matrix_banded sll_m_spline_matrix_banded module~sll_m_spline_matrix->module~sll_m_spline_matrix_banded module~sll_m_spline_matrix_base sll_m_spline_matrix_base module~sll_m_spline_matrix->module~sll_m_spline_matrix_base module~sll_m_spline_matrix_dense sll_m_spline_matrix_dense module~sll_m_spline_matrix->module~sll_m_spline_matrix_dense module~sll_m_spline_matrix->module~sll_m_working_precision module~sll_m_bsplines_base->module~sll_m_working_precision module~sll_m_bsplines_non_uniform->module~sll_m_bsplines_base module~sll_m_bsplines_non_uniform->module~sll_m_working_precision module~sll_m_bsplines_uniform->module~sll_m_bsplines_base module~sll_m_bsplines_uniform->module~sll_m_working_precision module~sll_m_spline_1d->module~sll_m_bsplines_base module~sll_m_spline_1d->module~sll_m_working_precision module~sll_m_spline_matrix_banded->iso_fortran_env module~sll_m_spline_matrix_banded->module~sll_m_spline_matrix_base module~sll_m_spline_matrix_banded->module~sll_m_working_precision module~sll_m_spline_matrix_base->module~sll_m_working_precision module~sll_m_spline_matrix_dense->iso_fortran_env module~sll_m_spline_matrix_dense->module~sll_m_spline_matrix_base module~sll_m_spline_matrix_dense->module~sll_m_working_precision

initialize the type sbase with polynomial degree, continuity ( -1: disc, 1: full) and number of gauss points per element

DO j=1,nD DO i=j+1,nD sf%A_axis(i,:,iBC)=sf%A_axis(i,:,iBC)-sf%A_axis(j,:,iBC)*(sf%A_axis(i,j,iBC)/sf%A_axis(j,j,iBC)) END DO END DO DO i=diri+1,nD sf%A_axis(i,:,iBC)=sf%A_axis(i,:,iBC)/sf%A_axis(i,i,iBC) END DO

NOT CLEAR HOW TO PRECONDITION A_EDGE WITH LU...

sf%R_Axis(1:nD,1:nD,iBC)=INV(sf%A_Axis(1:nD,1:nD,iBC)) sf%R_Axis(:,:,iBC)=TRANSPOSE(MATMUL(sf%R_Axis(:,:,iBC),sf%A_Axis(:,:,iBC))) DO i=1,deg+1; DO j=1,i-1 sf%R_Axis(i,j,iBC)=-sf%R_Axis(i,j,iBC) END DO; END DO sf%R_Edge(nBase-nD+1:nBase,nBase-nD+1:nBase,iBC)=INV(sf%A_Edge(nBase-nD+1:nBase,nBase-nD+1:nBase,iBC)) sf%R_Edge(:,:,iBC)=TRANSPOSE(MATMUL(sf%R_Edge(:,:,iBC),sf%A_Edge(:,:,iBC))) DO i=nBase-deg,nBase; DO j=i+1,nBase sf%R_Edge(i,j,iBC)=-sf%R_Edge(i,j,iBC) END DO; END DO

sf%R_axis( 1:nD ,:,iBC)=0.0_wp sf%R_edge(nBase-nD+1:nBase,:,iBC)=0.0_wp sf%R_Axis(1:nD,1:nD,iBC)=INV(sf%A_Axis(1:nD,1:nD,iBC)) sf%R_Axis( : , : ,iBC)=TRANSPOSE(MATMUL(sf%R_Axis(:,:,iBC),sf%A_Axis(:,:,iBC))) DO i=1,deg+1; DO j=1,i-1 sf%R_Axis(i,j,iBC)=-sf%R_Axis(i,j,iBC) END DO; END DO

sf%R_axis(1:nD,:,iBC)=0.0_wp

Type Bound

t_sBase

Arguments

Type IntentOptional Attributes Name
class(t_sBase), intent(inout) :: sf

self
integer, intent(in) :: deg_in

polynomial degree
integer, intent(in) :: continuity_in
class(t_sGrid), intent(in), TARGET :: grid_in

grid information
integer, intent(in) :: degGP_in

gauss quadrature points: nGP=degGP+1 per elements

Calls

proc~~sbase_init~~CallsGraph proc~sbase_init t_sBase%sBase_init add_element add_element proc~sbase_init->add_element barycentricweights barycentricweights proc~sbase_init->barycentricweights dmatip dmatip proc~sbase_init->dmatip eval_basis eval_basis proc~sbase_init->eval_basis eval_basis_and_n_derivs eval_basis_and_n_derivs proc~sbase_init->eval_basis_and_n_derivs eval_deriv eval_deriv proc~sbase_init->eval_deriv factorize factorize proc~sbase_init->factorize get_interp_points get_interp_points proc~sbase_init->get_interp_points init init proc~sbase_init->init initializevandermonde initializevandermonde proc~sbase_init->initializevandermonde legendregaussnodesandweights legendregaussnodesandweights proc~sbase_init->legendregaussnodesandweights mthpolynomialderivativematrix mthpolynomialderivativematrix proc~sbase_init->mthpolynomialderivativematrix proc~getlu getLU proc~sbase_init->proc~getlu proc~inv INV proc~sbase_init->proc~inv proc~sbase_alloc sBase_alloc proc~sbase_init->proc~sbase_alloc proc~sbase_test sBase_test proc~sbase_init->proc~sbase_test proc~sll_s_bsplines_new sll_s_bsplines_new proc~sbase_init->proc~sll_s_bsplines_new proc~sll_s_spline_matrix_new sll_s_spline_matrix_new proc~sbase_init->proc~sll_s_spline_matrix_new proc~solvemat SOLVEMAT proc~sbase_init->proc~solvemat swrite swrite proc~sbase_init->swrite xiip xiip proc~sbase_init->xiip dgetrf dgetrf proc~getlu->dgetrf proc~inv->dgetrf dgetri dgetri proc~inv->dgetri proc~sbase_alloc->dmatip proc~sbase_alloc->xiip wbaryip wbaryip proc~sbase_alloc->wbaryip proc~sbase_test->swrite proc~sbase_applybctodof_lgm t_sBase%sBase_applyBCtoDOF_LGM proc~sbase_test->proc~sbase_applybctodof_lgm proc~sbase_change_base t_sBase%sBase_change_base proc~sbase_test->proc~sbase_change_base proc~sbase_compare t_sBase%sBase_compare proc~sbase_test->proc~sbase_compare proc~sbase_eval t_sBase%sBase_eval proc~sbase_test->proc~sbase_eval proc~sbase_evaldof_base t_sBase%sBase_evalDOF_base proc~sbase_test->proc~sbase_evaldof_base proc~sbase_evaldof_gp t_sBase%sBase_evalDOF_GP proc~sbase_test->proc~sbase_evaldof_gp proc~sbase_evaldof_s t_sBase%sBase_evalDOF_s proc~sbase_test->proc~sbase_evaldof_s proc~sbase_initdof t_sBase%sBase_initDOF proc~sbase_test->proc~sbase_initdof proc~sbase_new sBase_new proc~sbase_test->proc~sbase_new proc~sll_s_bsplines_new->init sll_assert sll_assert proc~sll_s_bsplines_new->sll_assert proc~sll_s_spline_matrix_new->init sll_error sll_error proc~sll_s_spline_matrix_new->sll_error proc~solvemat->dgetrf dgetrs dgetrs proc~solvemat->dgetrs proc~solve SOLVE proc~sbase_applybctodof_lgm->proc~solve proc~sbase_change_base->proc~sbase_compare proc~sbase_change_base->proc~sbase_initdof eval eval proc~sbase_change_base->eval evalDOF_base evalDOF_base proc~sbase_change_base->evalDOF_base proc~sgrid_compare t_sGrid%sGrid_compare proc~sbase_compare->proc~sgrid_compare proc~sbase_eval->eval_basis proc~sbase_eval->eval_basis_and_n_derivs lagrangeinterpolationpolys lagrangeinterpolationpolys proc~sbase_eval->lagrangeinterpolationpolys proc~sgrid_find_elem t_sGrid%sGrid_find_elem proc~sbase_eval->proc~sgrid_find_elem proc~sbase_evaldof_s->proc~sbase_eval proc~sbase_evaldof_s->proc~sbase_evaldof_base compute_interpolant compute_interpolant proc~sbase_initdof->compute_interpolant proc~sbase_new->proc~sbase_init proc~solve->dgetrf proc~solve->dgetrs

Called by

proc~~sbase_init~~CalledByGraph proc~sbase_init t_sBase%sBase_init proc~sbase_test sBase_test proc~sbase_init->proc~sbase_test proc~sbase_copy t_sBase%sBase_copy proc~sbase_copy->proc~sbase_init proc~sbase_new sBase_new proc~sbase_new->proc~sbase_init proc~base_copy t_base%base_copy proc~base_copy->proc~sbase_copy proc~base_new Base_new proc~base_new->proc~sbase_new proc~readstatefilefromascii ReadStateFileFromASCII proc~readstatefilefromascii->proc~sbase_new proc~readstatefilefromascii->proc~base_new proc~sbase_test->proc~sbase_new interface~readstate ReadState interface~readstate->proc~readstatefilefromascii proc~init_base Init_Base proc~init_base->proc~base_new proc~initmhd3d t_functional_mhd3d%InitMHD3D proc~initmhd3d->proc~base_new proc~transform_sfl_init t_transform_sfl%transform_SFL_init proc~transform_sfl_init->proc~base_new proc~init_gvec_to_jorek init_gvec_to_jorek proc~init_gvec_to_jorek->interface~readstate proc~init_gvec_to_jorek->proc~init_base proc~restartfromstate RestartFromState proc~restartfromstate->interface~readstate proc~transform_sfl_new transform_sfl_new proc~transform_sfl_new->proc~transform_sfl_init

Source Code

SUBROUTINE sBase_init( sf,deg_in,continuity_in,grid_in,degGP_in)
! MODULES
USE MODgvec_GLobals,   ONLY: PI
USE MODgvec_LinAlg ,   ONLY: INV,SOLVEMAT,getLU
USE MODgvec_Basis1D,   ONLY:  LegendreGaussNodesAndWeights
USE MODgvec_Basis1D,   ONLY:  BarycentricWeights,InitializeVandermonde,MthPolynomialDerivativeMatrix
USE sll_m_bsplines,ONLY: sll_s_bsplines_new
USE sll_m_spline_interpolator_1d         ,ONLY: sll_t_spline_interpolator_1d
USE sll_m_boundary_condition_descriptors ,ONLY: sll_p_greville
USE sll_m_spline_matrix                  ,ONLY: sll_s_spline_matrix_new
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  INTEGER       , INTENT(IN   )        :: deg_in        !! polynomial degree
  INTEGER       , INTENT(IN   )        :: continuity_in !! continuity:
                                                        !! 0: disc. polynomial
                                                        !! deg-1: spline with cont. deg-1
  CLASS(t_sgrid), INTENT(IN   ),TARGET :: grid_in       !! grid information
  INTEGER       , INTENT(IN   )        :: degGP_in      !! gauss quadrature points: nGP=degGP+1 per elements
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  CLASS(t_sbase), INTENT(INOUT)        :: sf !! self
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  INTEGER  :: i,m,iGP,iElem,imin,jmin,nElems,deg,degGP,nBase
  INTEGER  :: iBC,j,diri,odd_even,nD,BC_derivs(deg_in+1),drv
  REAL(wp),ALLOCATABLE,DIMENSION(:,:) :: locbasis,VdmGP
  REAL(wp),DIMENSION(1:deg_in+1,1:deg_in+1) :: PLmat,Umat
!===================================================================================================================================
  IF(.NOT.test_called) THEN
    SWRITE(UNIT_stdOut,'(4X,A,3(A,I3),A)')'INIT sBase type:', &
         ' degree= ',deg_in, &
         ' gauss points per elem = ',degGP_in, &
         ' continuity= ',continuity_in, ' ...'
  END IF
  IF(sf%initialized) THEN
    CALL abort(__STAMP__, &
        "Trying to reinit sbase!")
  END IF
  IF(degGP_in.LT.deg_in) &
    CALL abort(__STAMP__, &
        "error in sbase: degGP must be > deg!")
  SELECT TYPE(sf)
  TYPE IS(t_sbase_disc)
    IF(continuity_in.NE.-1) &
      CALL abort(__STAMP__, &
          "error in sbase init: type is disc but continuity is not -1, mabye sbase_new was not called before!")
  TYPE IS(t_sbase_spl)
    IF(continuity_in.NE.deg_in-1) &
      CALL abort(__STAMP__, &
          "error in sbase init: type is spl but continuity is not deg-1, mabye sbase_new was not called before!")
  CLASS DEFAULT
      CALL abort(__STAMP__, &
          "error in sbase init: type is neither disc or spl!")
  END SELECT !Type
  sf%deg        =deg_in
  sf%continuity =continuity_in
  sf%grid       => grid_in
  sf%degGP      =  degGP_in

  nElems      = sf%grid%nElems
  deg         = sf%deg
  degGP       = sf%degGP

  sf%nGP  = (degGP+1)*nElems
  sf%nGP_str = (degGP+1)*(sf%grid%nElems_str-1)+1 !< for MPI
  sf%nGP_end = (degGP+1)*sf%grid%nElems_end

  IF(sf%continuity.EQ.-1)THEN !discontinuous
    sf%nBase  = (deg+1)*nElems
  ELSEIF(sf%continuity.EQ.deg-1)THEN !bspline with full continuity and interpolation base at boundaries
    sf%nBase  = nElems + deg
  ELSE
   CALL abort(__STAMP__, &
          'other spline continuities not yet implemented')
  END IF !continuity

  nBase = sf%nBase
  CALL sbase_alloc(sf)


  CALL LegendreGaussNodesAndWeights(degGP,sf%xi_GP,sf%w_GPloc) ![-1,1] !!!

  DO iElem=1,nElems
    sf%w_GP(1+(degGP+1)*(iElem-1):(degGP+1)*iElem)=                    0.5_wp *sf%w_GPloc(:)      *sf%grid%ds(iElem)
    sf%s_GP(1+(degGP+1)*(iElem-1):(degGP+1)*iElem)=sf%grid%sp(iElem-1)+0.5_wp*(sf%xi_GP(:)+1.0_wp)*sf%grid%ds(iElem)
  END DO !iElem

  SELECT TYPE(sf)
  TYPE IS(t_sbase_disc)
    IF(deg.EQ.0)THEN
      sf%xiIP(0)=0.0_wp
      CALL BarycentricWeights(deg,sf%xiIP,sf%wBaryIP)
      sf%DmatIP=0.0_wp
      DO iElem=1,nElems
        sf%base_offset(iElem)=1+(deg+1)*(iElem-1)
        sf%base_GP   (0:degGP,0,iElem)=1.0_wp
        sf%base_ds_GP(0:degGP,0,iElem)=0.0_wp
      END DO !iElem
      !zero deriv: evaluation of basis functions
      sf%base_dsAxis(0,1      )=1.0_wp
      sf%base_dsEdge(0,nBase  )=1.0_wp
    ELSE
      ALLOCATE(VdmGP( 0:degGP,0:deg))
      !  use chebychev-lobatto points for interpolation (closed form!), interval [-1,1]
      DO i=0,deg
        sf%xiIP(i)=-COS(REAL(i,wp)/REAL(deg,wp)*PI)
      END DO
      CALL BarycentricWeights(deg,sf%xiIP,sf%wBaryIP)
      CALL InitializeVandermonde(deg,degGP,sf%wBaryIP,sf%xiIP,sf%xi_GP,VdmGP)
      CALL MthPolynomialDerivativeMatrix(deg,sf%xiIP,1,sf%DmatIP)
      ! eval basis and  basis derivative
      DO iElem=1,nElems
        sf%base_offset(iElem)=1+(deg+1)*(iElem-1)
        sf%base_GP   (0:degGP,0:deg,iElem)=VdmGP(:,:)
        sf%base_ds_GP(0:degGP,0:deg,iElem)=MATMUL(VdmGP,sf%DmatIP)*(2.0_wp/sf%grid%ds(iElem))
      END DO !iElem
      !zero deriv: evaluation of basis functions (lagrange property!)
      sf%base_dsAxis(0,1      )=1.0_wp
      sf%base_dsAxis(0,2:deg+1)=0.0_wp
      sf%base_dsEdge(0,nBase-deg:nBase-1)=0.0_wp
      sf%base_dsEdge(0,          nBase  )=1.0_wp
      ! eval basis deriv at boundaries d/ds = d/dxi dxi/ds = 1/(0.5ds) d/dxi
      sf%base_dsAxis(1,1:deg+1        )=sf%DmatIP(  0,:)*(2.0_wp/sf%grid%ds(1))
      sf%base_dsEdge(1,nBase-deg:nBase)=sf%DmatIP(deg,:)*(2.0_wp/sf%grid%ds(nElems))
      !  higher derivatives
      DO i=2,deg
        sf%base_dsAxis(i,1:deg+1        )=MATMUL(TRANSPOSE(sf%DmatIP),sf%base_dsAxis(i-1,:))*(2.0_wp/sf%grid%ds(1))
        sf%base_dsEdge(i,nBase-deg:nBase)=MATMUL(TRANSPOSE(sf%DmatIP),sf%base_dsEdge(i-1,:))*(2.0_wp/sf%grid%ds(nElems))
      END DO
      DEALLOCATE(VdmGP)
    END IF !deg=0
    !interpolation:
    !  points are repeated at element interfaces (discontinuous)
    ALLOCATE(sf%s_IP(nBase)) !for spl, its allocated elsewhere...
    DO iElem=1,nElems
      sf%s_IP(1+(deg+1)*(iElem-1):(deg+1)*iElem)=sf%grid%sp(iElem-1)+0.5_wp*(sf%xiIP+1.0_wp)*sf%grid%ds(iElem)
    END DO !iElem
    sf%s_IP(1)=0.0_wp
    sf%s_IP(nBase)=1.0_wp
  TYPE IS(t_sbase_spl)
    ALLOCATE(locbasis(0:deg,0:deg))
    CALL sll_s_bsplines_new(sf%bspl ,degree=deg,periodic=.FALSE., &
                            xmin=sf%grid%sp(0),xmax=sf%grid%sp(nElems),&
                            ncells=nElems,breaks=sf%grid%sp)
    !basis evaluation

    IF(sf%bspl%nBasis.NE.nBase) CALL abort(__STAMP__, &
                                           'problem with bspl basis')
    DO iElem=1,nElems
      j=1+(degGP+1)*(iElem-1)
      CALL sf%bspl % eval_basis(sf%s_GP(j),sf%base_GP(0,0:deg,iElem),imin)
      IF(imin.EQ.-1) CALL abort(__STAMP__,&
                                'problem, element in eval_basis not found!')
      DO iGP=1,degGP
        CALL sf%bspl % eval_basis(sf%s_GP(j+iGP),sf%base_GP(iGP,0:deg,iElem),jmin)
        IF(jmin.EQ.-1) CALL abort(__STAMP__,&
                                 'problem, element in eval_basis not found!')
        IF(jmin.NE.imin) CALL abort(__STAMP__,&
                                    'problem, GP are not in one element!')
      END DO !iGP=0,degGP
      CALL sf%bspl % eval_deriv(sf%s_GP(j),sf%base_ds_GP(0,0:deg,iElem),imin)
      IF(imin.EQ.-1) CALL abort(__STAMP__,&
                                'problem, element in eval_basis not found!')
      DO iGP=1,degGP
        CALL sf%bspl % eval_deriv(sf%s_GP(j+iGP),sf%base_ds_GP(iGP,0:deg,iElem),jmin)
        IF(jmin.EQ.-1) CALL abort(__STAMP__,&
                                  'problem, element in eval_basis not found!')
        IF(jmin.NE.imin) CALL abort(__STAMP__,&
                                    'problem, GP are not in one element!')
      END DO !iGP=0,degGP
      sf%base_offset(iElem)=imin
    END DO !iElem=1,nElems
    !eval all basis derivatives at boundaries
    CALL sf%bspl % eval_basis_and_n_derivs(sf%grid%sp(     0),deg,locBasis,imin) !locBasis(0:nderiv,0:deg Base)
    IF(imin.NE.1) CALL abort(__STAMP__,&
                             'problem eval_deriv left')
    sf%base_dsAxis(0:deg,1:deg+1) =locbasis(:,:) ! basis functions 1 ...deg+1

    CALL sf%bspl % eval_basis_and_n_derivs(sf%grid%sp(nElems),deg,locbasis,imin)
    IF(imin.NE.nBase-deg) CALL abort(__STAMP__,&
                                     'problem eval_deriv right')
    sf%base_dsEdge(0:deg,nBase-deg:nBase)=locbasis(:,:) ! basis functions nBase-deg ... nbase

    !interpolation
    CALL sf%Interpol%init (sf%bspl,sll_p_greville,sll_p_greville)
    CALL sf%Interpol%get_interp_points ( sf%s_IP )

    DEALLOCATE(locbasis)
  END SELECT !TYPE is t_sbase_disc/spl


  !mass matrix
  CALL sll_s_spline_matrix_new(sf%mass , "banded",nBase,sf%deg,sf%deg)
  DO iElem=1,nElems
    jmin=sf%base_offset(iElem)
    DO i=0,deg
      DO j=0,deg
        CALL sf%mass%add_element(jmin+i,jmin+j, &
             (SUM(sf%w_GP(1+(degGP+1)*(iElem-1):(degGP+1)*iElem)*sf%base_GP(:,i,iElem)*sf%base_GP(:,j,iElem))))
      END DO !j=0,deg
    END DO !i=0,deg
  END DO !iElem=1,nElems
  CALL sf%mass % factorize()


  !STRONG BOUNDARY CONDITIONS: A and R matrices

  DO iBC=1,NBC_TYPES
    SELECT CASE(iBC)      !nDOF involved:    Dirichlet?  odd(0),even(1)
    CASE(BC_TYPE_OPEN)     ; nD =0                                   ! do nothing
    CASE(BC_TYPE_NEUMANN)  ; nD =1          ;   diri=0 ;  odd_even=0 ! first derivative=0
    CASE(BC_TYPE_DIRICHLET); nD =1          ;   diri=1 ;  odd_even=0 ! odd_even not used, first DOF=BC_val
    CASE(BC_TYPE_SYMM)     ; nD =(deg+1)/2  ;   diri=0 ;  odd_even=0 ! open,    all derivatives (2*k-1)=0, k=1,...(deg+1)/2
    CASE(BC_TYPE_SYMMZERO) ; nD =1+(deg+1)/2;   diri=1 ;  odd_even=0 ! dirichlet=0+ derivatives (2*k-1)=0, k=1,...(deg+1)/2
    CASE(BC_TYPE_ANTISYMM) ; nD =1+deg/2    ;   diri=1 ;  odd_even=1 ! dirichlet=0+ derivatives (2*k  )=0  k=1,... deg/2
    END SELECT !iBC
    sf%nDOF_BC(iBC)=nD
    !A and R are already initialized as unit matrices!!
    IF((nD.GT.0).AND.(deg.GT.0))THEN
      IF(diri.EQ.1) BC_derivs(1)=0  !A(1,:) is already unit matrix
      DO i=diri+1,nD
         BC_derivs(i)=2*(i-diri)-(1-odd_even)
      END DO
      DO i=diri+1,nD
        drv=BC_derivs(i)
        sf%A_Axis(        i,:,iBC)=sf%base_dsAxis(drv,:)/sf%base_dsAxis(drv,i) !normalized with diagonal entry
        sf%A_Edge(nBase+1-i,:,iBC)=sf%base_dsEdge(drv,:)/sf%base_dsEdge(drv,nBase+1-i)
      END DO

      !PRECONDITIONING OF AMAT, FORWARD SWEEP:

      !lower left triangular part to zero, and normalize diagonal
      !! DO j=1,nD
      !!   DO i=j+1,nD
      !!     sf%A_axis(i,:,iBC)=sf%A_axis(i,:,iBC)-sf%A_axis(j,:,iBC)*(sf%A_axis(i,j,iBC)/sf%A_axis(j,j,iBC))
      !!   END DO
      !! END DO
      !! DO i=diri+1,nD
      !!   sf%A_axis(i,:,iBC)=sf%A_axis(i,:,iBC)/sf%A_axis(i,i,iBC)
      !! END DO

      ! A*x=0 , [A_s A_r]*[x1 x2] = 0  A_s*x1=-A_r*x2, U^T (PL)^T x1 = -A_r*x2 -> NEW A_s= (PL)^T, A_r = U^{-T} A_r
      CALL getLU(nD,TRANSPOSE(sf%A_axis(1:nD,1:nD,iBC)),PLmat(1:nD,1:nD),Umat(1:nD,1:nD))
      sf%A_axis(1:nD,1:nD      ,iBC)=TRANSPOSE(PLmat(1:nD,1:nD))
      sf%A_axis(1:nD,nD+1:deg+1,iBC)=SOLVEMAT(TRANSPOSE(Umat(1:nD,1:nD)),sf%A_axis(1:nD,nD+1:deg+1,iBC))

      !upper right triangular part to zero, and normalize diagonal
      DO j=1,nD
        DO i=j+1,nD
          sf%A_Edge(nBase+1-i,:,iBC)= sf%A_Edge(nBase+1-i,:,iBC)  &
                                     -sf%A_Edge(nBase+1-j,:,iBC)*(sf%A_Edge(nBase+1-i,nBase+1-j,iBC)/sf%A_Edge(nBase+1-j,nBase+1-j,iBC))
        END DO
      END DO
      DO i=diri+1,nD
        sf%A_Edge(nBase+1-i,:,iBC)=sf%A_Edge(nBase+1-i,:,iBC)/sf%A_Edge(nBase+1-i,nBase+1-i,iBC)
      END DO

      !! NOT CLEAR HOW TO PRECONDITION A_EDGE WITH LU...



      !invert BC part (A_s is ( nD x nD ) and invertible):
      !     |  A_s   A_r |        |  0                    0 |
      !  A= |            |  , R=  |                         |
      !     |   0     I  |        | -(A_s^{-1}A_r)^T     I  |
      !
      sf%R_axis(   1:nD   ,:   ,iBC)=0.0_wp
      sf%R_axis(nD+1:deg+1,1:nD,iBC)=-TRANSPOSE(SOLVEMAT(sf%A_axis(1:nD,1:nD,iBC),sf%A_axis(1:nD,nD+1:deg+1,iBC)))

      sf%AR_axis(      1:nD,:,iBC)=sf%A_axis(      1:nD,:,iBC)
      sf%AR_axis(nD+1:deg+1,:,iBC)=sf%R_axis(nD+1:deg+1,:,iBC)
      !     |   I     0  |        |  I    -(B_s^{-1}B_r)^T  |
      !  B= |            |  , R=  |                         |
      !     |  B_r   B_s |        |  0                    0 |
      !
      sf%R_edge(nBase-deg:nBase-nD,nBase+1-nD:nBase,iBC)=-TRANSPOSE(SOLVEMAT(sf%A_edge(nBase+1-nD:nBase,nBase+1-nD:nBase  ,iBC), &
                                                                             sf%A_edge(nBase+1-nD:nBase,nBase-deg:nBase-nD,iBC)))
      sf%R_edge(nBase+1-nD:nBase,:,iBC)=0.0_wp


      sf%AR_edge(nBase-deg:nBase-nD,:,iBC)= sf%R_edge(nBase-deg:nBase-nD,:,iBC)
      sf%AR_edge(nBase+1-nD:nBase  ,:,iBC)= sf%A_edge(nBase+1-nD:nBase  ,:,iBC)


      !! sf%R_Axis(1:nD,1:nD,iBC)=INV(sf%A_Axis(1:nD,1:nD,iBC))
      !! sf%R_Axis(:,:,iBC)=TRANSPOSE(MATMUL(sf%R_Axis(:,:,iBC),sf%A_Axis(:,:,iBC)))
      !! DO i=1,deg+1; DO j=1,i-1
      !!   sf%R_Axis(i,j,iBC)=-sf%R_Axis(i,j,iBC)
      !! END DO; END DO
      !! sf%R_Edge(nBase-nD+1:nBase,nBase-nD+1:nBase,iBC)=INV(sf%A_Edge(nBase-nD+1:nBase,nBase-nD+1:nBase,iBC))
      !! sf%R_Edge(:,:,iBC)=TRANSPOSE(MATMUL(sf%R_Edge(:,:,iBC),sf%A_Edge(:,:,iBC)))
      !! DO i=nBase-deg,nBase; DO j=i+1,nBase
      !!   sf%R_Edge(i,j,iBC)=-sf%R_Edge(i,j,iBC)
      !! END DO; END DO
      !prepare for applyBC
      sf%invA_axis(:,:,iBC)=INV(sf%A_axis(:,:,iBC))
      sf%invA_edge(:,:,iBC)=INV(sf%A_edge(:,:,iBC))
      !automatically set rows 1:nD to zero for R matrices (no contribution from these DOF)
      !! sf%R_axis(         1:nD   ,:,iBC)=0.0_wp
      !! sf%R_edge(nBase-nD+1:nBase,:,iBC)=0.0_wp
    END IF
  END DO!iBC=1,NBC_TYPES

  ! mode number dependent boundary condition at axis, m  <=> -iBC
  ! iBC= 0, m=0:  BC_TYPE_SYMM
  ! iBC=-1, m=1:  BC_TYPE_ANTISYMM
  ! iBC=-2, m=2:  BC_TYPE_SYMMZERO
  ! ... odd m<=deg: ANTISYMM + all derivatives 1,...,m-1 =0
  ! ...even m<=deg: SYMMZERO + all derivatives 1,...,m-1 =0
  !m=0
  sf%nDOF_BC(   0) = sf%nDOF_BC(   BC_TYPE_SYMM)
  sf%A_Axis(:,:,0) = sf%A_Axis(:,:,BC_TYPE_SYMM)

  !  odd m <=deg
  DO m=1,deg+1,2
    iBC=-m
    !BC_TYPE_ANTISYMM: nD =1+deg/2    ;   diri=1 ;  odd_even=1 ! dirichlet=0+ derivatives (2*k  )=0  k=1,... deg/2
    !and also derivatives (2*k-1)=0, k=1,...,(m-1)/2  (2*k-1 < m-1)
    nD=1+deg/2 +(m-1)/2
    sf%nDOF_BC(iBC)=nD
    IF(nD.EQ.deg+1) CYCLE ! All deg+1 DOFs at the boundary just set to zero (A_axis is already identity matrix), if deg even and m=deg+1
    i=1
    BC_derivs(1)=0 !diri
    DO drv=1,(m-1)
      i=i+1
      BC_derivs(i)=drv
    END DO
    DO drv=m+1,deg,2
      i=i+1
      BC_derivs(i)=drv
    END DO
    IF(i.NE.nD) CALL abort(__STAMP__,&
                           "DEBUG odd m, i NE nD")
    DO i=2,nD
      drv=BC_derivs(i)
      sf%A_Axis(i,:,iBC)=sf%base_dsAxis(drv,:)/sf%base_dsAxis(drv,i) !normalized with diagonal entry
    END DO
  END DO
  !  even m <=deg
  DO m=2,deg+1,2
    iBC=-m
    !BC_TYPE_SYMMZERO: nD =1+(deg+1)/2;   diri=1 ;  odd_even=0 ! dirichlet=0+ derivatives (2*k-1)=0, k=1,...(deg+1)/2
    !and also derivatives (2*k)=0, k=1,...,(m-2)/2  (2*k < m-1)
    nD=1+(deg+1)/2 +(m-2)/2
    sf%nDOF_BC(iBC)=nD
    IF(nD.EQ.deg+1) CYCLE ! All deg+1 DOFs at the boundary just set to zero (A_axis is already identity matrix), if deg odd and m=deg+1
    i=1
    BC_derivs(1)=0 !diri
    DO drv=1,(m-1)
      i=i+1
      BC_derivs(i)=drv
    END DO
    DO drv=m+1,deg,2
      i=i+1
      BC_derivs(i)=drv
    END DO
    IF(i.NE.nD) CALL abort(__STAMP__,&
                           "DEBUG even m,i NE nD")
    DO i=2,nD
      drv=BC_derivs(i)
      sf%A_Axis(i,:,iBC)=sf%base_dsAxis(drv,:)/sf%base_dsAxis(drv,i) !normalized with diagonal entry
    END DO
  END DO

  DO m=0,deg+1
    iBC=-m
    nD=sf%nDOF_BC(iBC)
    IF(nD.LT.deg+1)THEN
      !precondition A via LU decomp. of A_s^T
      ! A*x=0 , [A_s A_r]*[x1 x2] = 0  A_s*x1=-A_r*x2, U^T (PL)^T x1 = -A_r*x2 -> NEW A_s= (PL)^T, A_r = U^{-T} A_r
      CALL getLU(nD,TRANSPOSE(sf%A_axis(1:nD,1:nD,iBC)),PLmat(1:nD,1:nD),Umat(1:nD,1:nD))
      sf%A_axis(1:nD,1:nD      ,iBC)=TRANSPOSE(PLmat(1:nD,1:nD))
      sf%A_axis(1:nD,nD+1:deg+1,iBC)=SOLVEMAT(TRANSPOSE(Umat(1:nD,1:nD)),sf%A_axis(1:nD,nD+1:deg+1,iBC))

      !invert BC part (A_s is ( nD x nD ) and invertible):
      !     |  A_s   A_r |        |  0                    0 |
      !  A= |            |  , R=  |                         |
      !     |   0     I  |        | -(A_s^{-1}A_r)^T     I  |
      !
      sf%R_axis(   1:nD   ,:   ,iBC)=0.0_wp
      sf%R_axis(nD+1:deg+1,1:nD,iBC)=-TRANSPOSE(SOLVEMAT(sf%A_axis(1:nD,1:nD,iBC),sf%A_axis(1:nD,nD+1:deg+1,iBC)))

      sf%AR_axis(      1:nD,:,iBC)=sf%A_axis(      1:nD,:,iBC)
      sf%AR_axis(nD+1:deg+1,:,iBC)=sf%R_axis(nD+1:deg+1,:,iBC)
      !invert BC part
      !!sf%R_Axis(1:nD,1:nD,iBC)=INV(sf%A_Axis(1:nD,1:nD,iBC))
      !!sf%R_Axis( :  , :  ,iBC)=TRANSPOSE(MATMUL(sf%R_Axis(:,:,iBC),sf%A_Axis(:,:,iBC)))
      !!DO i=1,deg+1; DO j=1,i-1
      !!  sf%R_Axis(i,j,iBC)=-sf%R_Axis(i,j,iBC)
      !!END DO; END DO
      !prepare for applyBC
      sf%invA_axis(:,:,iBC)=INV(sf%A_axis(:,:,iBC))
      !automatically set rows 1:nD to zero for R matrices (no contribution from these DOF)
      !!sf%R_axis(1:nD,:,iBC)=0.0_wp
    ELSE
      !nD=deg+1, A_axis is simply the identity!
      sf%R_axis(   :,:,iBC)=0.0_wp
      sf%AR_axis(  :,:,iBC)=sf%A_axis(:,:,iBC)
      sf%invA_axis(:,:,iBC)=sf%A_axis(:,:,iBC)
    END IF
  END DO!m=0,deg+1


  sf%initialized=.TRUE.
  IF(.NOT.test_called) THEN
    SWRITE(UNIT_stdOut,'(4X,A)')'... DONE'
    CALL sBase_test(sf)
  END IF

END SUBROUTINE sBase_init