test sbase variable
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(t_sBase), | intent(inout) | :: | sf |
self |
SUBROUTINE sBase_test( sf) ! MODULES USE MODgvec_GLobals, ONLY: testdbg,testlevel,nfailedMsg,nTestCalled,testUnit IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES CLASS(t_sBase), INTENT(INOUT) :: sf !! self !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES INTEGER :: i,iTest,iElem,jElem,BC_Type(2),deg,degGP,nElems,cont,nBase REAL(wp) :: x,y,y2,dy,dy2 REAL(wp) :: y_BC(0:sf%deg),y2_BC(0:sf%deg),base_x(0:sf%deg) REAL(wp) :: g_IP(1:sf%nBase),dofs(1:sf%nBase) REAL(wp) :: y_GP(1:sf%nGP) REAL(wp) :: dof_GP(1:sf%nGP),BC_Val(2) REAL(wp),PARAMETER :: realtol=1.0E-11_wp CHARACTER(LEN=10) :: fail CLASS(t_sbase),ALLOCATABLE :: testsbase LOGICAL :: check(5) REAL(wp),ALLOCATABLE :: oldDOF(:,:),newDOF(:,:) !=================================================================================================================================== test_called=.TRUE. IF(testlevel.LE.0) RETURN IF(testdbg) THEN Fail=" DEBUG !!" ELSE Fail=" FAILED !!" END IF deg=sf%deg degGP=sf%degGP cont=sf%continuity nBase=sf%nbase nElems=sf%grid%nElems nTestCalled=nTestCalled+1 SWRITE(UNIT_stdOut,'(A,I4,A)')'>>>>>>>>> RUN SBASE TEST ID',nTestCalled,' >>>>>>>>>' IF(testlevel.GE.1)THEN iTest=101 ; IF(testdbg)WRITE(*,*)'iTest=',iTest IF(testdbg.OR.(.NOT.( ABS(sf%s_IP(1)).LT. realtol))) THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(2(A,I4),(A,E11.3))') & ' degree = ',deg, ' continuity = ',cont, & '\n => should be 0.0 : sIP(1)= ',sf%s_IP(1) END IF !TEST iTest=102 ; IF(testdbg)WRITE(*,*)'iTest=',iTest IF(testdbg.OR.(.NOT.( ABS(sf%s_IP(nBase)-1.0_wp).LT. realtol))) THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),(A,E11.3))') & ' degree = ',deg, & ' continuity = ',cont, ', nBase= ',nBase, & '\n => should be 1.0 : sIP(nBase)= ',sf%s_IP(nBase) END IF !TEST iTest=103 ; IF(testdbg)WRITE(*,*)'iTest=',iTest IF(testdbg.OR.(.NOT.( ABS(SUM(sf%w_GP)-1.0_wp).LT. realtol))) THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),(A,E11.3))') & ' degree = ',deg, & ' continuity = ',cont, ', nBase= ',nBase, & '\n => should be 1.0 : SUM(w_GP)= ',SUM(sf%w_GP) END IF !TEST iTest=104 ; IF(testdbg)WRITE(*,*)'iTest=',iTest CALL sf%eval(0.0_wp,0,iElem,base_x) IF(testdbg.OR.(.NOT.((ABS(base_x(0)-1.).LT.realtol).AND. & ( SUM(ABS(base_x(1:deg))).LT. realtol).AND. & (iElem.EQ.1) ))) THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),A,I6,*(E11.3))') & ' degree = ',deg, & ' continuity = ',cont, ', nBase= ',nBase, & '\n => should be 1 and (1,0,0...) : base_x(x=0)= ',iElem,base_x END IF !TEST iTest=105 ; IF(testdbg)WRITE(*,*)'iTest=',iTest CALL sf%eval(1.0_wp,0,iElem,base_x) IF(testdbg.OR.(.NOT.((ABS(base_x(deg)-1.) .LT. realtol ).AND. & (MAXVAL(ABS(base_x(0:deg-1))).LT. realtol ).AND. & (iElem .EQ. nElems ) ))) THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),A,*(E11.3))') & ' degree = ',deg, & ' continuity = ',cont, ', nBase= ',nBase, & '\n => should be ...,0,0,1 : base_x(x=1)= ',base_x END IF !TEST iTest=106 ; IF(testdbg)WRITE(*,*)'iTest=',iTest jElem=(nElems+1)/2 x=0.5_wp*(sf%grid%sp(jElem-1)+sf%grid%sp(jElem)) CALL sf%eval(x,0,iElem,base_x) IF(testdbg.OR.(.NOT.(iElem.EQ.jElem)))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),2(A,I6))') & ' degree = ',deg, & ' continuity = ',cont, ', nBase= ',nBase, & '\n => should be ',jElem,': iElem= ' , iElem END IF !TEST iTest=107 ; IF(testdbg)WRITE(*,*)'iTest=',iTest jElem=MIN(2,nElems) x=sf%grid%sp(jElem-1)+0.01_wp*sf%grid%ds(jElem) CALL sf%eval(x ,0,iElem,base_x) IF(testdbg.OR.(.NOT.(iElem.EQ.jElem)))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),2(A,I6))') & ' degree = ',deg, & ' continuity = ',cont, ', nBase= ',nBase, & '\n => should be ',jElem,': iElem= ' , iElem END IF !TEST !check compare iTest=111 ; IF(testdbg)WRITE(*,*)'iTest=',iTest CALL sbase_new(testsbase,sf%deg,sf%continuity,sf%grid, sf%degGP) CALL testsbase%compare(sf,is_same=check(1)) IF(.NOT.check(1))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),(A))') & ' degree = ',deg, & ' continuity = ',cont, ', nBase= ',nBase, & '\n => should be true' END IF !TEST !check compare iTest=112 ; IF(testdbg)WRITE(*,*)'iTest=',iTest CALL testsbase%compare(sf, cond_out=check(1:5)) IF(.NOT.ALL(check(:)))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),(A,5L))') & ' degree = ',deg, & ' continuity = ',cont, ', nBase= ',nBase, & '\n => should be all true', check(:) END IF !TEST CALL testsbase%free() DEALLOCATE(testsbase) !check compare iTest=113 ; IF(testdbg)WRITE(*,*)'iTest=',iTest CALL sbase_new(testsbase,sf%deg+1,MERGE(-1,sf%deg,(sf%continuity.EQ.-1)),sf%grid, sf%degGP) CALL testsbase%compare(sf,is_same=check(1)) IF(check(1))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),(A))') & ' degree = ',deg, & ' continuity = ',cont, ', nBase= ',nBase, & '\n => should be false' END IF !TEST CALL testsbase%free() DEALLOCATE(testsbase) !check change_base, execution (aborts if not correct) iTest=121 ; IF(testdbg)WRITE(*,*)'iTest=',iTest CALL sbase_new(testsbase,sf%deg,sf%continuity,sf%grid, sf%degGP) ALLOCATE(oldDOF(1:sf%nBase,2),newDOF(1:testsBase%nBase,2)) oldDOF(:,1)=1.0_wp oldDOF(:,2)=2.0_wp CALL testsbase%change_base(sf, 2, oldDOF,newDOF) y = MAXVAL(ABS(newDOF(:,1)-oldDOF(:,1))) y2 = MINVAL(ABS(newDOF(:,2)-oldDOF(:,2))) IF(testdbg.OR.(.NOT.((y .LT.realtol ).AND. & (y2 .LT.realtol ))))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),2(A,E11.3))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,'\n => should be zero ...',y & ,'\n => should be zero ...',y2 END IF !TEST DEALLOCATE(oldDOF,newDOF) CALL testsbase%free() DEALLOCATE(testsbase) !check change_base,only execution (aborts if not correct) iTest=122 ; IF(testdbg)WRITE(*,*)'iTest=',iTest CALL sbase_new(testsbase,sf%deg+1,MERGE(-1,sf%deg,(sf%continuity.EQ.-1)),sf%grid, sf%degGP) ALLOCATE(oldDOF(2,1:sf%nBase),newDOF(2,1:testsBase%nBase)) oldDOF(1,:)=-1.0_wp oldDOF(2,:)=-2.0_wp CALL testsbase%change_base(sf, 1, oldDOF,newDOF) y = MINVAL(ABS(newDOF(1,:)))-MINVAL(ABS(oldDOF(1,:))) y2= MINVAL(ABS(newDOF(2,:)))-MINVAL(ABS(oldDOF(2,:))) DEALLOCATE(oldDOF,newDOF) CALL testsbase%free() DEALLOCATE(testsbase) IF(testdbg.OR.(.NOT.((y .LT.realtol ).AND. & (y2 .LT.realtol ))))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),2(A,E11.3))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,'\n => should be zero ...',y & ,'\n => should be zero ...',y2 END IF !TEST END IF !testlevel>=1 IF(testlevel.GE.2)THEN jElem=MAX(1,nElems/2) x=sf%grid%sp(jElem-1) IF((cont.EQ.-1).AND.(nElems.GT.1))THEN g_IP = testf(sf%s_IP(:)) g_IP(1:(deg+1)*jElem) =g_IP(1:(deg+1)*jElem)-0.113_wp !discontinuous between jElem and jElem+1 ELSE g_IP = testf(sf%s_IP) END IF !TEST dofs(1:sf%nBase)=sf%initDOF(g_IP) iTest=201 ; IF(testdbg)WRITE(*,*)'iTest=',iTest !CHECK boundary values y=testf(0.0_wp) IF((cont.EQ.-1).AND.(nElems.GT.1))THEN y=y-0.113_wp END IF y2=testf(1.0_wp) IF(testdbg.OR.(.NOT.((ABS( dofs(1)-y ).LT.realtol ).AND.& (ABS( dofs(nBase)-y2 ).LT.realtol ))))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),4(A,E11.3))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,'\n => should be ',y ,' : dofs(1) = ' , dofs(1) & ,'\n => should be ',y2 ,' : dofs(nBase) = ' , dofs(nBase) END IF !TEST !check dsAxis and dsEdge, basis evaluation iTest=202 ; IF(testdbg)WRITE(*,*)'iTest=',iTest y=SUM(sf%base_dsAxis(0,:)*dofs(1:deg+1)) IF((cont.EQ.-1).AND.(nElems.GT.1))THEN y=y+0.113_wp END IF y2=SUM(sf%base_dsEdge(0,:)*dofs(nbase-deg:nBase)) IF(testdbg.OR.(.NOT.((ABS( y - testf(0.0_wp) ).LT.realtol ).AND.& (ABS( y2 - testf(1.0_wp) ).LT.realtol ))))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),4(A,E11.3))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,'\n => should be ',testf(0.0_wp) ,' : dofs(y=0) = ' , y & ,'\n => should be ',testf(1.0_wp) ,' : dofs(y=1) = ' , y2 END IF !TEST !check dsAxis and dsEdge, first derivative iTest=203 ; IF(testdbg)WRITE(*,*)'iTest=',iTest y =sf%evalDOF_s(0.0_wp,1,dofs) y2=sf%evalDOF_s(1.0_wp,1,dofs) IF(deg.GE.1)THEN dy =SUM(sf%base_dsAxis(1,:)*dofs( 1:deg+1)) dy2=SUM(sf%base_dsEdge(1,:)*dofs(nbase-deg:nBase)) IF(testdbg.OR.(.NOT.((ABS( dy - testf_dx(0.0_wp) ).LT.realtol/sf%grid%ds( 1) ).AND.& (ABS( dy2 - testf_dx(1.0_wp) ).LT.realtol/sf%grid%ds(nElems) ).AND.& (ABS( y - dy ).LT.realtol/sf%grid%ds( 1) ).AND.& (ABS( y2 - dy2 ).LT.realtol/sf%grid%ds(nElems) ))))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),8(A,E11.3))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,'\n => should be ',testf_dx(0.0_wp) ,' : dofs_ds(y=0) = ' , dy & ,'\n => should be ',testf_dx(1.0_wp) ,' : dofs_ds(y=1) = ' , dy2 & ,'\n => should be ',dy ,' : dofs_ds(y=0) = ' , y & ,'\n => should be ',dy2 ,' : dofs_ds(y=1) = ' , y2 END IF !TEST ELSE IF(testdbg.OR.(.NOT.((ABS( y ).LT.realtol ).AND.& (ABS( y2 ).LT.realtol ))))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),4(A,E11.3))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,'\n => should be ',0.0_wp ,' : dofs_ds^2(y=0) = ' , dy & ,'\n => should be ',0.0_wp ,' : dofs_ds^2(y=1) = ' , dy2 END IF !TEST END IF !check dsAxis and dsEdge, second derivative iTest=204 ; IF(testdbg)WRITE(*,*)'iTest=',iTest dy =sf%evalDOF_s(0.0_wp,2,dofs) dy2=sf%evalDOF_s(1.0_wp,2,dofs) IF(deg.GE.2)THEN y =SUM(sf%base_dsAxis(2,:)*dofs(1:deg+1)) y2=SUM(sf%base_dsEdge(2,:)*dofs(nbase-deg:nBase)) IF(testdbg.OR.(.NOT.((ABS( dy - testf_dxdx(0.0_wp) ).LT.realtol/(sf%grid%ds( 1)**2) ).AND.& (ABS( dy2 - testf_dxdx(1.0_wp) ).LT.realtol/(sf%grid%ds(nElems)**2) ).AND.& (ABS( y - dy ).LT.realtol/(sf%grid%ds( 1)**2) ).AND.& (ABS( y2 - dy2 ).LT.realtol/(sf%grid%ds(nElems)**2) ))))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),8(A,E11.3))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,'\n => should be ',testf_dxdx(0.0_wp) ,' : dofs_ds^2(y=0) = ' , dy & ,'\n => should be ',testf_dxdx(1.0_wp) ,' : dofs_ds^2(y=1) = ' , dy2 & ,'\n => should be ',dy ,' : dofs_ds(y=0) = ' , y & ,'\n => should be ',dy2 ,' : dofs_ds(y=1) = ' , y2 END IF !TEST ELSE IF(testdbg.OR.(.NOT.((ABS( dy ).LT.realtol ).AND.& (ABS( dy2 ).LT.realtol ))))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),4(A,E11.3))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,'\n => should be ',0.0_wp ,' : dofs_ds^2(y=0) = ' , dy & ,'\n => should be ',0.0_wp ,' : dofs_ds^2(y=1) = ' , dy2 END IF !TEST END IF IF(nElems.GE.2)THEN !check eval, evalDOF_base and evalDOF_s, function evaluation iTest=211 ; IF(testdbg)WRITE(*,*)'iTest=',iTest x=sf%grid%sp(jElem)+0.9503_wp*sf%grid%ds(jElem+1) !element jElem+1 CALL sf%eval(x ,0,iElem,base_x) y=sf%evalDOF_base(iElem,base_x,dofs(:)) y2=sf%evalDOF_s(x,0,dofs(:)) IF(testdbg.OR.(.NOT.((ABS( y-testf(x) ).LT.realtol ).AND.& (ABS( y-y2 ).LT.realtol ).AND.& (iElem .EQ.jElem+1 ) )))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),2(A,I6),4(A,E11.3))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,'\n => should be ',jElem,': iElem= ' , iElem & ,'\n => should be ',testf(x),' : y = ' , y & ,'\n => should be ',y ,' : y2= ' , y2 END IF !TEST END IF ! nElems.GE.2 !check eval, evalDOF_base and evalDOF_s, function evaluation, in jElem iTest=212 ; IF(testdbg)WRITE(*,*)'iTest=',iTest x=sf%grid%sp(jElem-1)+ 0.7353_wp*sf%grid%ds(jElem) !in element jElem CALL sf%eval(x ,0,iElem,base_x) y=sf%evalDOF_base(iElem,base_x,dofs(:)) y2=sf%evalDOF_s(x,0,dofs(:)) IF((cont.EQ.-1).AND.(nElems.GT.1))THEN y=y+0.113_wp y2=y2+0.113_wp END IF IF(testdbg.OR.(.NOT.((ABS( y-testf(x) ).LT.realtol ).AND.& (ABS( y-y2 ).LT.realtol ).AND.& (iElem .EQ.jElem ) )))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),2(A,I6),4(A,E11.3))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,'\n => should be ',jElem,': iElem= ' , iElem & ,'\n => should be ',testf(x),' : y = ' , y & ,'\n => should be ',y ,' : y2= ' , y2 END IF !test IF(nElems.GE.2)THEN !test first derivative iTest=213 ; IF(testdbg)WRITE(*,*)'iTest=',iTest x=sf%grid%sp(jElem)+0.64303_wp*sf%grid%ds(jElem+1) !in elem jElem+1 CALL sf%eval(x ,1,iElem,base_x) dy=sf%evalDOF_base(iElem,base_x,dofs(:)) dy2=sf%evalDOF_s(x,1,dofs(:)) IF(testdbg.OR.(.NOT.((ABS(dy-testf_dx(x)).LT.realtol/sf%grid%ds(jElem+1) ).AND.& (ABS(dy-dy2 ).LT.realtol/sf%grid%ds(jElem+1) ).AND.& (iElem .EQ.jElem+1 ))))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),2(A,I6),4(A,E11.3))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,'\n => should be ',jElem,': iElem= ' , iElem & ,'\n => should be ',testf_dx(x),': dy = ' , dy & ,'\n => should be ',y ,': dy2= ' , dy END IF !test !test second derivative iTest=214 ; IF(testdbg)WRITE(*,*)'iTest=',iTest x=sf%grid%sp(jElem)+0.17313_wp*sf%grid%ds(jElem+1) !in elem jElem+1 CALL sf%eval(x ,2,iElem,base_x) dy=sf%evalDOF_base(iElem,base_x,dofs(:)) !second derivative dy2=sf%evalDOF_s(x,2,dofs(:)) IF(testdbg.OR.(.NOT.((ABS(dy-testf_dxdx(x)).LT.realtol/(sf%grid%ds(jElem+1)**2) ).AND.& (ABS(dy-dy2 ).LT.realtol/(sf%grid%ds(jElem+1)**2) ).AND.& (iElem .EQ.jElem+1 ))))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),2(A,I6),4(A,E11.3))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,'\n => should be ',jElem,': iElem= ' , iElem & ,'\n => should be ',testf_dxdx(x),': dy = ' , dy & ,'\n => should be ',y ,': dy2= ' , dy2 END IF !test END IF ! nElems.GE.2 !check evalDOF_GP iTest=221 ; IF(testdbg)WRITE(*,*)'iTest=',iTest y_GP = testf(sf%s_GP) IF((cont.EQ.-1).AND.(nElems.GT.1)) y_GP(1:(degGP+1)*jElem)=y_GP(1:(degGP+1)*jElem)-0.113_wp dof_GP = sf%evalDOF_GP(0,dofs) IF(testdbg.OR.(.NOT.((MAXVAL(ABS(y_GP-dof_GP)).LT.realtol ))))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),4(A,E11.3))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,'\n => should be ',y_GP(1) ,': dof_GP(1) = ' , dof_GP(1) & ,'\n => should be ',MAXVAL(ABS(y_GP)) ,': MAX(|dof_GP|) = ' , MAXVAL(ABS(dof_GP)) END IF !test !check integration, full domain iTest=222 ; IF(testdbg)WRITE(*,*)'iTest=',iTest y = testf_int(0.0_wp,1.0_wp) IF((cont.EQ.-1).AND.(nElems.GT.1)) y=y-0.113_wp*(sf%grid%sp(jElem)-sf%grid%sp(0)) y2 = SUM(dof_GP*sf%w_GP) IF(testdbg.OR.(.NOT.( (ABS(y-y2).LT.realtol ) )))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),2(A,E11.3))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,'\n => should be ',y ,': SUM(dof_GP*wGP) = ' , y2 END IF !test !check integration, last element iTest=223 ; IF(testdbg)WRITE(*,*)'iTest=',iTest y = testf_int(sf%grid%sp(nElems-1),1.0_wp) y2 = SUM( dof_GP(1+(degGP+1)*(nElems-1):(degGP+1)*nElems) & *sf%w_GP(1+(degGP+1)*(nElems-1):(degGP+1)*nElems)) IF(testdbg.OR.(.NOT.( (ABS(y-y2).LT.realtol ) )))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),2(A,E11.3))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,'\n => should be ',y ,': SUM(dof_GP*wGP)|_lastElem = ' , y2 END IF !test !check evalDOF_GP, first derivative iTest=224 ; IF(testdbg)WRITE(*,*)'iTest=',iTest y_GP = testf_dx(sf%s_GP) dof_GP = sf%evalDOF_GP(DERIV_S,dofs) IF(testdbg.OR.(.NOT.((MAXVAL(ABS(y_GP-dof_GP)).LT.realtol/MINVAL(sf%grid%ds) ))))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),4(A,E11.3))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,'\n => should be ',y_GP(1) ,': dof_ds_GP(1) = ' , dof_GP(1) & ,'\n => should be ',MAXVAL(ABS(y_GP)) ,': MAX(|dof_ds_GP|) = ' , MAXVAL(ABS(dof_GP)) END IF !test !apply BC iTest=231 ; IF(testdbg)WRITE(*,*)'iTest=',iTest g_IP = dofs ! BC_Type(1)=BC_TYPE_DIRICHLET BC_Type(2)=BC_TYPE_OPEN BC_Val(1:2)=(/0.96_wp,0.63_wp/) CALL sf%applyBCtoDOF(g_IP,BC_Type,BC_val) y_BC=g_IP(1:deg+1) y2_BC( 0)=0.96_wp y2_BC(1:deg)=dofs(2:deg+1) IF(testdbg.OR.(.NOT.((MAXVAL(ABS(y_BC-y2_BC)).LT.realtol ))))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),4(A,E11.3))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,'\n => should be ',y_BC(0) ,': y(BCaxis) = ' , y2_BC(0) & ,'\n => should be ',MAXVAL(ABS(y_BC)) ,': MAX(|y(BC_axis)|) = ' , MAXVAL(ABS(y2_BC)) END IF !test iTest=232 ; IF(testdbg)WRITE(*,*)'iTest=',iTest g_IP = dofs ! BC_Type(1)=BC_TYPE_OPEN BC_Type(2)=BC_TYPE_DIRICHLET CALL sf%applyBCtoDOF(g_IP,BC_type,BC_Val) y_BC=g_IP(nBase-deg:nBase) y2_BC(0:deg-1)=dofs(nBase-deg:nBase-1) y2_BC( deg)=0.63_wp IF(testdbg.OR.(.NOT.((MAXVAL(ABS(y_BC-y2_BC)).LT.realtol ))))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),4(A,E11.3))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,'\n => should be ',y_BC(deg) ,': y(BCedge) = ' , y2_BC(deg) & ,'\n => should be ',MAXVAL(ABS(y_BC)) ,': MAX(|y(BCedge)|) = ' , MAXVAL(ABS(y2_BC)) END IF !test iTest=233 ; IF(testdbg)WRITE(*,*)'iTest=',iTest g_IP = dofs ! BC_Type(1)=BC_TYPE_NEUMANN BC_Type(2)=BC_TYPE_OPEN BC_Val(1:2)=(/0.96_wp,0.63_wp/) CALL sf%applyBCtoDOF(g_IP,BC_Type,BC_val) y_BC=MATMUL(sf%base_dsAxis(:,:),g_IP(1:deg+1)) DO i=1,deg y_BC(i)=y_BC(i)*sf%grid%ds(1)**i/REAL(i*i,wp) END DO IF(testdbg.OR.(.NOT.((ABS(y_BC(1)).LT.realtol ))))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),A,*(E13.5))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,"\n => should be y,y',y''... : dy/ds BC_NEUMANN axis = " , y_BC(:) END IF !test iTest=234 ; IF(testdbg)WRITE(*,*)'iTest=',iTest g_IP = dofs ! BC_Type(1)=BC_TYPE_OPEN BC_Type(2)=BC_TYPE_NEUMANN CALL sf%applyBCtoDOF(g_IP,BC_type,BC_Val) y_BC=MATMUL(sf%base_dsEdge(:,:),g_IP(nBase-deg:nBase)) DO i=1,deg y_BC(i)=y_BC(i)*sf%grid%ds(nElems)**i/REAL(i*i,wp) END DO IF(testdbg.OR.(.NOT.((ABS(y_BC(1)).LT.realtol ))))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),A,*(E13.5))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,"\n => should be y,y',y'',... : dy/ds BC_NEUMANN edge = " , y_BC(:) END IF !test iTest=235 ; IF(testdbg)WRITE(*,*)'iTest=',iTest g_IP = dofs ! BC_Type(1)=BC_TYPE_SYMM BC_Type(2)=BC_TYPE_OPEN BC_Val(1:2)=(/0.96_wp,0.63_wp/) CALL sf%applyBCtoDOF(g_IP,BC_Type,BC_val) y_BC=MATMUL(sf%base_dsAxis(:,:),g_IP(1:deg+1)) DO i=1,deg y_BC(i)=y_BC(i)*sf%grid%ds(1)**i/REAL(i*i,wp) END DO IF(testdbg.OR.(.NOT.((MAXVAL(ABS(y_BC(1:deg:2))).LT.realtol ))))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),A,*(E13.5))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,"\n => should be y,y',y''... : dy/ds BC_SYMM axis = " , y_BC(:) END IF !test iTest=236 ; IF(testdbg)WRITE(*,*)'iTest=',iTest g_IP = dofs ! BC_Type(1)=BC_TYPE_OPEN BC_Type(2)=BC_TYPE_SYMM CALL sf%applyBCtoDOF(g_IP,BC_type,BC_Val) y_BC=MATMUL(sf%base_dsEdge(:,:),g_IP(nBase-deg:nBase)) DO i=1,deg y_BC(i)=y_BC(i)*sf%grid%ds(nElems)**i/REAL(i*i,wp) END DO IF(testdbg.OR.(.NOT.((MAXVAL(ABS(y_BC(1:deg:2))).LT.realtol ))))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),A,*(E13.5))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,"\n => should be y,y',y'',... : dy/ds BC_SYMM edge = " , y_BC(:) END IF !test iTest=237 ; IF(testdbg)WRITE(*,*)'iTest=',iTest g_IP = dofs ! BC_Type(1)=BC_TYPE_SYMMZERO BC_Type(2)=BC_TYPE_OPEN BC_Val(1:2)=(/0.96_wp,0.63_wp/) CALL sf%applyBCtoDOF(g_IP,BC_Type,BC_val) y_BC=MATMUL(sf%base_dsAxis(:,:),g_IP(1:deg+1)) DO i=1,deg y_BC(i)=y_BC(i)*sf%grid%ds(1)**i/REAL(i*i,wp) END DO IF(testdbg.OR.(.NOT.((MAXVAL(ABS(y_BC(1:deg:2))).LT.realtol ).AND. & (ABS(y_BC(0) ).LT.realtol ) )))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),A,*(E13.5))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,"\n => should be y,y',y''... : dy/ds BC_SYMMZERO axis = " , y_BC(:) END IF !test iTest=238 ; IF(testdbg)WRITE(*,*)'iTest=',iTest g_IP = dofs ! BC_Type(1)=BC_TYPE_OPEN BC_Type(2)=BC_TYPE_SYMMZERO CALL sf%applyBCtoDOF(g_IP,BC_type,BC_Val) y_BC=MATMUL(sf%base_dsEdge(:,:),g_IP(nBase-deg:nBase)) DO i=1,deg y_BC(i)=y_BC(i)*sf%grid%ds(nElems)**i/REAL(i*i,wp) END DO IF(testdbg.OR.(.NOT.((MAXVAL(ABS(y_BC(1:deg:2))).LT.realtol ).AND. & (ABS(y_BC(0) ).LT.realtol ) )))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),A,*(E13.5))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,"\n => should be y,y',y'',... : dy/ds BC_SYMMZERO edge = " , y_BC(:) END IF !test iTest=239 ; IF(testdbg)WRITE(*,*)'iTest=',iTest g_IP = dofs ! BC_Type(1)=BC_TYPE_ANTISYMM BC_Type(2)=BC_TYPE_OPEN BC_Val(1:2)=(/0.96_wp,0.63_wp/) CALL sf%applyBCtoDOF(g_IP,BC_Type,BC_val) y_BC=MATMUL(sf%base_dsAxis(:,:),g_IP(1:deg+1)) DO i=1,deg y_BC(i)=y_BC(i)*(sf%grid%ds(1)**i)/REAL(i*i,wp) END DO IF(testdbg.OR.(.NOT.((MAXVAL(ABS(y_BC(2:deg:2))).LT.realtol ).AND. & (ABS(y_BC(0) ).LT.realtol ) )))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),A,*(E13.5))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,"\n => should be y,y',y''... : dy/ds BC_ANTISYMM axis = " , y_BC(:) END IF !test iTest=240 ; IF(testdbg)WRITE(*,*)'iTest=',iTest g_IP = dofs ! BC_Type(1)=BC_TYPE_OPEN BC_Type(2)=BC_TYPE_ANTISYMM CALL sf%applyBCtoDOF(g_IP,BC_type,BC_Val) y_BC=MATMUL(sf%base_dsEdge(:,:),g_IP(nBase-deg:nBase)) DO i=1,deg y_BC(i)=y_BC(i)*sf%grid%ds(nElems)**i/REAL(i*i,wp) END DO IF(testdbg.OR.(.NOT.((MAXVAL(ABS(y_BC(2:deg:2))).LT.realtol ).AND. & (ABS(y_BC(0) ).LT.realtol ) )))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),A,*(E13.5))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,"\n => should be y,y',y'',... : dy/ds BC_ANTISYMM edge = " , y_BC(:) END IF !test iTest=241 ; IF(testdbg)WRITE(*,*)'iTest=',iTest g_IP = dofs ! BC_Type(1)=-deg ! TEST m-dependent BC_AXIS...iBC=-m, m=deg <=> function and all derivatives up to deg-1 are zero BC_Type(2)=BC_TYPE_OPEN CALL sf%applyBCtoDOF(g_IP,BC_type,BC_Val) y_BC=MATMUL(sf%base_dsAxis(:,:),g_IP(1:deg+1)) DO i=1,deg y_BC(i)=y_BC(i)*(sf%grid%ds(1)**i)/REAL(i*i,wp) END DO IF(testdbg.OR.(.NOT.((MAXVAL(ABS(y_BC(1:deg-1))).LT.realtol ).AND. & (ABS(y_BC(0) ).LT.realtol ) )))THEN nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') & '\n!! SBASE TEST ID',nTestCalled ,': TEST ',iTest,Fail nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(3(A,I4),A,*(E13.5))') & ' degree = ',deg & ,' continuity = ',cont, ', nBase= ',nBase & ,"\n => should be y, dy^k/ds, k=1,deg-1 at axis = " , y_BC(:) END IF !test END IF !testlevel>=2 test_called=.FALSE. CONTAINS !FUNCTION WITH CONTINUITY deg-1 at sp(jElem) (discontinuity must be treated on an array level, see above) ELEMENTAL FUNCTION testf(x) REAL(wp),INTENT(IN) :: x REAL(wp) :: testf IF(sf%deg.EQ.0) THEN IF(x.GT.sf%grid%sp(jElem))THEN testf=0.21_wp ELSE testf=0.05_wp END IF ELSE testf=(1.13_wp-0.3_wp*x)**sf%deg & + MAX(0.0_wp,0.21_wp*(x-sf%grid%sp(jElem)))**sf%deg END IF END FUNCTION testf ELEMENTAL FUNCTION testf_dx(x) REAL(wp),INTENT(IN) :: x REAL(wp) :: testf_dx IF(sf%deg.EQ.0)THEN testf_dx=0.0_wp ELSEIF(sf%deg.EQ.1)THEN IF(x.GT.sf%grid%sp(jElem))THEN testf_dx=-0.3_wp*REAL(sf%deg,wp)+ 0.21_wp*REAL(sf%deg,wp) ELSE testf_dx=-0.3_wp*REAL(sf%deg,wp) END IF ELSE testf_dx=-0.3_wp*REAL(sf%deg,wp)*(1.13_wp-0.3_wp*x)**(sf%deg-1) & + 0.21_wp*REAL(sf%deg,wp)*MAX(0.0_wp,0.21_wp*(x-sf%grid%sp(jElem)))**(sf%deg-1) END IF END FUNCTION testf_dx ELEMENTAL FUNCTION testf_dxdx(x) REAL(wp),INTENT(IN) :: x REAL(wp) :: testf_dxdx IF(sf%deg.LE.1)THEN testf_dxdx=0.0_wp ELSEIF(sf%deg.EQ.2)THEN IF(x.GT.sf%grid%sp(jElem))THEN testf_dxdx= (0.3_wp)**2*REAL(sf%deg*(sf%deg-1),wp) & + (0.21_wp)**2*REAL(sf%deg*(sf%deg-1),wp) ELSE testf_dxdx=(0.3_wp)**2*REAL(sf%deg*(sf%deg-1),wp) END IF ELSE testf_dxdx=(0.3_wp)**2*REAL(sf%deg*(sf%deg-1),wp)*(1.13_wp-0.3_wp*x)**(sf%deg-2) & + (0.21_wp)**2*REAL(sf%deg*(sf%deg-1),wp)*MAX(0.0_wp,0.21_wp*(x-sf%grid%sp(jElem)))**(sf%deg-2) END IF END FUNCTION testf_dxdx FUNCTION testf_int(a,b) REAL(wp),INTENT(IN) :: a,b !!a<=b REAL(wp) :: testf_int IF(a.GT.sf%grid%sp(jElem))THEN IF(sf%deg.EQ.0)THEN testf_int=0.21_wp*(b-a) ELSE testf_int= 1.0_wp/(-0.3_wp*REAL(sf%deg+1,wp))*( (1.13_wp-0.3_wp*b)**(sf%deg+1) & -(1.13_wp-0.3_wp*a)**(sf%deg+1)) & + 1.0_wp/(0.21_wp*REAL(sf%deg+1,wp))*( (0.21_wp*(b-sf%grid%sp(jElem)))**(sf%deg+1) & -(0.21_wp*(a-sf%grid%sp(jElem)))**(sf%deg+1) ) END IF ELSEIF(b.GE.sf%grid%sp(jElem))THEN IF(sf%deg.EQ.0)THEN testf_int=0.21_wp*(b-sf%grid%sp(jElem))+0.05_wp*(sf%grid%sp(jElem)-a) ELSE testf_int= 1.0_wp/(-0.3_wp*REAL(sf%deg+1,wp))*( (1.13_wp-0.3_wp*b)**(sf%deg+1) & -(1.13_wp-0.3_wp*a)**(sf%deg+1)) & + 1.0_wp/(0.21_wp*REAL(sf%deg+1,wp))*( (0.21_wp*(b-sf%grid%sp(jElem)))**(sf%deg+1) ) END IF ELSE IF(sf%deg.EQ.0)THEN testf_int=0.05_wp*(b-a) ELSE testf_int= 1.0_wp/(-0.3_wp*REAL(sf%deg+1,wp))*( (1.13_wp-0.3_wp*b)**(sf%deg+1) & -(1.13_wp-0.3_wp*a)**(sf%deg+1)) END IF END IF END FUNCTION testf_int END SUBROUTINE sbase_test