C USER INPUT FOR ADAPTIVE MESH CONSTRAINT C SUBROUTINE UMESHMOTION(UREF,ULOCAL,NODE,NNDOF, $ LNODETYPE,ALOCAL,NDIM,TIME,DTIME,PNEWDT, $ KSTEP,KINC,KMESHSWEEP,JMATYP,JGVBLOCK,LSMOOTH) C include 'ABA_PARAM.INC' C C USER DEFINED DIMENSION STATEMENTS C CHARACTER*80 PARTNAME C The dimensions of the variables ARRAY C must be set equal to or greater than 15 DIMENSION ARRAY(1000),JPOS(15) DIMENSION ULOCAL(*) DIMENSION UGLOBAL(NDIM) DIMENSION JGVBLOCK(*),JMATYP(*) DIMENSION ALOCAL(NDIM,*) DIMENSION WVLOCAL(3),WVGLOBAL(3) PARAMETER (NELEMMAX=100) PARAMETER (CHARLENGTH = 1.D-1,ELINC = 0.05D0) PARAMETER (ASMALL = 1.D-10) DIMENSION JELEMLIST(NELEMMAX),JELEMTYPE(NELEMMAX) C C Wear topology common block and parameters C C nstreamlines = number of points at the reference section C nGenElem = number of general sectors in the model C nCylElem = number of cylindrical sectors in the model C (Cylindrical elements are not supported currently) C nRevOffset = node offset specified under SMG,revolve C nReflOffset = node offset specified under SMG,reflect C (Set it to zero if the model is not reflected) C jslnodes = The first component of this arry is the node number C at the reference section. The second component is the C Node thatprovides the wear direction. If the wear has C to be applied in the normal direction, set it to zero. C endisp1 = Array of cumulative ablation magnitude for streamlines C from the previous increment C endisp0 = Array of Cumulative ablation magnitude for streamlines C from the current increment C enstrlen1 = Array of streamline lenghts from the previous increment C enstrlen0 = Array of streamline lengths from the current increment C parameter (nstreamlines=27,nGenElem=48,nCylElem=0,nRevOffset=1000, $ nReflOffset=50000) common /wear/ $ jslnodes(2,nstreamlines), $ endisp0(nstreamlines),endisp1(nstreamlines), $ enstrlen0(nstreamlines),enstrlen1(nstreamlines), $ lvalidinc,lvalidsweep data lvalidinc /-1/ data jslnodes/ $40, 0, $41, 0, $42, 0, $11, 0, $43, 0, $44, 0, $45, 0, $46, 0, $12,47, $19,63, $64, 0, $65, 0, $66, 0, $20, 0, $50040, 0, $50041, 0, $50042, 0, $50011, 0, $50043, 0, $50044, 0, $50045, 0, $50046, 0, $50012,50047, $50019,50063, $50064, 0, $50065, 0, $50066, 0/ C NELEMS = NELEMMAX ldebug = 0 jtyp = 0 CALL GETNODETOELEMCONN(NODE,NELEMS,JELEMLIST,JELEMTYPE, $ JRCD,JGVBLOCK) luseendisp = 1 C if (ldebug.ne.0) then C write (7,*) 'UMESHMOTION' C write (7,*) 'JELEMLIST ' C do k1 = 1,NELEMS C write (7,*) JELEMLIST(k1) C end do C write (7,*) 'kstep ',kstep C write (7,*) 'kinc ',kinc C write (7,*) 'kmeshsweep ',kmeshsweep C write (7,*) 'nstreamlines ',nstreamlines C write (7,*) 'uref ',uref C C end if IF (KINC.EQ.1.AND.kmeshsweep.eq.0.and.lvalidinc.eq.-1) then C first time in this routine altogether lvalidinc = 2 lvalidsweep = 1 luseendisp = 0 else if (kinc .eq. lvalidinc) then C first time in this routine this new increment C Move back copy of energy to front copy do k1 = 1,nstreamlines endisp1(k1) = endisp0(k1) endisp0(k1) = 0.0d0 enstrlen1(k1)= enstrlen0(k1) enstrlen0(k1)= 0.0d0 end do lvalidinc = kinc + 1 lvalidsweep = 1 end if if (kmeshsweep.eq.lvalidsweep) then C first time in this routine this mesh sweep C reset back copy of energy dissipation and streamline length do k1 = 1,nstreamlines endisp0(k1) = 0.0d0 enstrlen0(k1) = 0.0d0 enstrlen1(k1)= enstrlen0(k1) enstrlen0(k1)= 0.0d0 end do lvalidsweep = kmeshsweep + 1 end if if (ldebug.ne.0) then write (7,*) 'lvalidinc ',lvalidinc write (7,*) 'lvalidsweep ',lvalidsweep end if LOCNUM = 0 JRCD = 0 PARTNAME = ' ' PEEQ = 0.0D0 CALL GETPARTINFO(NODE,0,PARTNAME,LOCNUM,JRCD) lstreamline = 0 lposition = 0 nslnodes=nGenElem+2*nCylElem C Which streamline is locnum on? do ksl = 1,nstreamlines C locnum > nReflOffset, node lies on the reflected half if ((locnum-nReflOffset).gt.0) then C process only jslnodes that are higher than nReflOffset if ((jslnodes(1,ksl)-nReflOffset).gt.0) then if (MOD((locnum-jslnodes(1,ksl)),nRevOffset).eq.0) then lstreamline=ksl lposition=1+(locnum-jslnodes(1,ksl))/nRevOffset end if end if else C locnum < nReflOffset, process only jslnodes less than nReflOffset if ((jslnodes(1,ksl)-nReflOffset).lt.0) then if (MOD((locnum-jslnodes(1,ksl)),nRevOffset).eq.0) then lstreamline=ksl lposition=1+(locnum-jslnodes(1,ksl))/nRevOffset end if end if end if end do if (ldebug.ne.0) then write (7,*) locnum, node, ' is on streamline ',lstreamline write (7,*) locnum, node, ' is at position ',lposition end if if (lstreamline.eq.0) return CALL GETVRMAVGATNODE(LOCNUM,JTYP,'CSTRESS',ARRAY,JRCD, $ JELEMLIST,NELEMS,JMATYP,JGVBLOCK) CPRESS = ARRAY(1) CSHEAR = SQRT(ARRAY(2)**2+ARRAY(3)**2) C IF (CSHEAR.LT.ASMALL) GOTO 500 CALL GETVRMAVGATNODE(LOCNUM,JTYP,'CDISP',ARRAY,JRCD, $ JELEMLIST,NELEMS,JMATYP,JGVBLOCK) CSLIP = SQRT(ARRAY(2)**2+ARRAY(3)**2) C IF (CSLIP.LT.ASMALL) GOTO 500 C Length of the streamline segment is 0.5*(distance_to_previous_node+distance_to_next_node) C C coordinates of locnum LTRN=0 CALL GETVRN(LOCNUM,'COORD',ARRAY,JRCD,JGVBLOCK,LTRN) current_x=ARRAY(1) current_y=ARRAY(2) current_z=ARRAY(3) C C coordinates of the previous point LPREV=0 if (lposition.eq.1) then LPREV=LOCNUM+(nslnodes-1)*nRevOffset else LPREV=LOCNUM-nRevOffset end if CALL GETVRN(LPREV,'COORD',ARRAY,JRCD,JGVBLOCK,LTRN) prev_x=ARRAY(1) prev_y=ARRAY(2) prev_z=ARRAY(3) C C coordinates of the next point... LNEXT=0 if (lposition.eq.nslnodes) then LNEXT=LOCNUM-(nslnodes-1)*nRevOffset else LNEXT=LOCNUM+nRevOffset end if CALL GETVRN(LNEXT,'COORD',ARRAY,JRCD,JGVBLOCK,LTRN) cnext_x=ARRAY(1) cnext_y=ARRAY(2) cnext_z=ARRAY(3) C C distances dist_prev= SQRT ((prev_x-current_x)**2+(prev_y-current_y)**2 $ +(prev_z-current_z)**2) dist_next=SQRT((cnext_x-current_x)**2+(cnext_y-current_y)**2+ $ (cnext_z-current_z)**2) SL_LENGTH=0.5*(dist_prev+dist_next) c c This must also include the length of the streamline segment if (ldebug.ne.0) then write (7,*) 'Adding ',(SL_LENGTH), $ ' to streamline ',lstreamline write (7,*) 'Adding ',(CSLIP * CSHEAR * $ SL_LENGTH),' to streamline ',lstreamline end if C ENSTRLEN0(lstreamline) = ENSTRLEN0(lstreamline) + $ (SL_LENGTH) C ENDISP0(lstreamline) = ENDISP0(lstreamline) + $ (CSLIP * CSHEAR * SL_LENGTH) C 500 if (luseendisp.ne.0) then ldebug2 = 0 if (ldebug.ne.0)ldebug2 = 1 if (ldebug.ne.0) then write (7,*) 'Ablating node ', $ locnum, ENDISP1(lstreamline) * UREF write (7,*) 'alocal' write (7,100) alocal(1,1),alocal(1,2),alocal(1,3) write (7,100) alocal(2,1),alocal(2,2),alocal(2,3) write (7,100) alocal(3,1),alocal(3,2),alocal(3,3) 100 format (1x,3(1pg15.2)) end if C if (ABS(enstrlen1(lstreamline)).lt.ASMALL) then SURFV=0.0d0 else SURFV = ENDISP1(lstreamline) * UREF / enstrlen1(lstreamline) end if C C Wear Direction: C Ablate the node along ULOCAL(3) if there is no dependency. C Else, compute the wear direction in the global system, apply the wear C and transform it back to the local system C if (jslnodes(2,lstreamline).eq.0) then ULOCAL(NDIM)=ULOCAL(NDIM)-SURFV else C Compute the node number of the wear master Master=(lposition-1)*nRevOffset+jslnodes(2,lstreamline) LTRN=0 CALL GETVRN(MASTER,'COORD',ARRAY,JRCD,JGVBLOCK,LTRN) cmaster_x=ARRAY(1) cmaster_y=ARRAY(2) cmaster_z=ARRAY(3) dist_master=SQRT((cmaster_x-current_x)**2+ $ (cmaster_y-current_y)**2+(cmaster_z-current_z)**2) C Wear in Global directions WVGLOBAL(1)=SURFV*(cmaster_x-current_x)/dist_master WVGLOBAL(2)=SURFV*(cmaster_y-current_y)/dist_master WVGLOBAL(3)=SURFV*(cmaster_z-current_z)/dist_master C Transformation to the local directions do k2=1,NDIM WVLOCAL(k2)=0 do k3=1,NDIM WVLOCAL(k2)=WVLOCAL(k2)+WVGLOBAL(k3)*ALOCAL(k3,k2) end do end do do k2 = 1,ndim ULOCAL(k2) = ULOCAL(k2) + WVLOCAL(k2) end do end if C C SHARE OF ELEMENT CONSUMED IN NEXT INCREMENT PELEM = DTIME * SURFV / CHARLENGTH PNEWDT0 = PNEWDT PNEWDT1 = 0.0D0 IF (PELEM.GT.ELINC.AND.SURFV.GT.0.0D0) THEN PNEWDT1 = CHARLENGTH * ELINC / SURFV IF (PNEWDT1.LT.PNEWDT0) THEN WRITE (7,*) 'CHANGING TIME INCREMENT FROM ',PNEWDT0 PNEWDT = PNEWDT1 WRITE (7,*) 'TO ',PNEWDT WRITE (7,*) 'BASED ON NODE ',LOCNUM END IF END IF end if C RETURN END