1! 2! CalculiX - A 3-dimensional finite element program 3! Copyright (C) 1998-2021 Guido Dhondt 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU General Public License as 7! published by the Free Software Foundation(version 2); 8! 9! 10! This program is distributed in the hope that it will be useful, 11! but WITHOUT ANY WARRANTY; without even the implied warranty of 12! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13! GNU General Public License for more details. 14! 15! You should have received a copy of the GNU General Public License 16! along with this program; if not, write to the Free Software 17! Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 18! 19! MASSLESS DYNAMIC CONTACT 20 21C> *MASSLESS DYNAMIC CONTACT*: Extracting the submatrices (i/b) from the global stiffness matrix 22 23C> Submatrices: 24C> + the bb matrix (neqtot x neqtot) 25C> + the bi matrix (neqtot x neq(1)) 26C> + the ib matrix (neq(1) x neqtot) 27C> + the ii matrix (neq(1) x neq(1)) 28C> - this matrix has the same 29C> size of the global matrix, in which the bi-entries, 30C> ib-entries, bb-entries (off-diagonal) are set to zero 31C> the diagonal bb-entries are set to one 32 33C> @param au LOWER triangle of STIFFNESS matrix of size: neq[0] 34C> @param ad DIAGONAL of STIFFNESS matrix of size: neq[0] 35C> @param jq Location in field **irow** of the first subdiagonal nonzero in column i (only for symmetric matrices) 36C> @param irow Row of element i in field au (i.e. au(i)) 37C> @param neq NTOT of equations: 38C> + neq[0]: mechanical TOTDOF, 39C> + neq[1]: TOTDOF_mech + TOTDOF_thermal 40C> + neq[2]: neq[1] + # of single point constraints (only for modal calculations) 41C> @param nzs see **massless.c** 42C> @param aubb the param 43C> @param adbb the param 44C> @param jqbb the param 45C> @param irowbb the param 46C> @param neqtot the param 47C> @param nzsbb the param 48C> @param islavnode see **massless.c** 49C> @param nslavs see **massless.c** 50C> @param imastnode see **massless.c** 51C> @param nmasts NTOT of master nodes 52C> @param aubi the param 53C> @param jqbi the param 54C> @param irowbi the param 55C> @param nzsbi the param 56C> @param auib the param 57C> @param jqib the param 58C> @param irowib the param 59C> @param nzsib the param 60C> @param kslav slave contact dofs list 61C> @param ktot slave and master dofs list 62 63C> @see massless.c 64C> @author Guido Dhondt 65 66 subroutine extract_matrices(au,ad,jq,irow,neq,nzs,aubb,adbb,jqbb, 67 & irowbb,neqtot,nzsbb,islavnode,nslavs,imastnode,nmasts,aubi, 68 & jqbi,irowbi,nzsbi,auib,jqib,irowib,nzsib,kslav,ktot, icolbb) 69 implicit none 70! 71 integer jq(*),irow(*),neq(*),nzs(*),jqbb(*),irowbb(*),neqtot, 72 & nzsbb,islavnode(*),nslavs,imastnode(*),nmasts,jqbi(*), 73 & irowbi(*),nzsbi,jqib(*),irowib(*),nzsib,kbb,i,j,kslav(*), 74 & ktot(*),id,lbb, icolbb(*) 75! 76 real*8 au(*),ad(*),aubb(*),adbb(*),aubi(*),auib(*) 77! 78! number of potential nonzero subdiagonal entries 79! 80 nzsbb=0 81 nzsbi=0 82 nzsib=0 83! 84! treating column by column in the global stiffness matrix 85! 86 kbb=1 87 do i=1,neq(1) 88C write(*,*)'DEBUG i',i 89! 90! filling submatrix bi 91! 92 if(i.ne.ktot(kbb)) then !does not correspond to MAS or SLV node therefore potentially to bi? 93 jqbi(i)=nzsbi+1 94 95! write(*,*)'DEBUG jq(i+1)-1',jq(i+1)-1 96 do j=jq(i),jq(i+1)-1 ! all items in column i in big K 97! write(*,*)'DEBUG j',j 98 call nident(ktot,irow(j),neqtot,id) ! if belongs to ktot (then is b) or not (then is i) entry 99 if(id.gt.0) then 100 if(ktot(id).eq.irow(j)) then ! ths its bi entry 101! write(*,*)'(its bi) i,j,ktot(id),: ',i,j,ktot(id) 102 nzsbi = nzsbi+1 103 aubi(nzsbi) = au(j) 104 au(j) = 0.d0 105 irowbi(nzsbi) = id 106 endif 107 endif 108 enddo 109 cycle ! was in bi, therefore we cycle 110 endif 111! then its a b column 112! filling submatrix ib and bb 113! 114 lbb=1 115 jqbi(i)=nzsbi+1! check if workssssss 116 jqbb(kbb)=nzsbb+1 117 jqib(kbb)=nzsib+1 ! 118! 119! LOOP over all NONDIAG entries of current columni of i of K matrix 120 loop: do j=jq(i),jq(i+1)-1 !QUESTION why the -1 ? 121! write(*,*)'DEBUG j',j 122 if(irow(j).lt.ktot(lbb)) then ! its an i row, not a b row. SMALLER 123! write(*,*)'(its ib (<) )', irow(j),ktot(lbb) 124 nzsib=nzsib+1 125 auib(nzsib)=au(j) 126 au(j)=0.d0 127 irowib(nzsib)=irow(j) ! TODO: off by +18 or more 128 cycle 129 elseif(irow(j).eq.ktot(lbb)) then ! b row its only if its equal EQUAL 130! write(*,*)'(its bb (=))', irow(j),ktot(lbb) 131 nzsbb=nzsbb+1 132 aubb(nzsbb)=au(j) 133 au(j)=0.d0 134 irowbb(nzsbb)=lbb 135 lbb=lbb+1 136 cycle 137 else ! BIGGER ( prob wont happen but just in case) --> happens all the time D: 138! write(*,*) 'bigger >! ',irow(j),ktot(lbb) 139 do !QUESTION is this a loop? 140 lbb=lbb+1 141! write(*,*) 'lbb++: ' , lbb ! debug 142 if(irow(j).lt.ktot(lbb)) then 143! write(*,*)'(its ib (<))', irow(j),ktot(lbb) 144 nzsib=nzsib+1 145 auib(nzsib)=au(j) 146 au(j)=0.d0 ! QUESTION why??? 147 irowib(nzsib)=irow(j) ! TODO: off by +18 or more 148 cycle loop ! TODO not sure if w/wo loop ?? 149 elseif(irow(j).eq.ktot(lbb)) then 150! write(*,*)'(its bb (=))', irow(j),ktot(lbb) 151 nzsbb=nzsbb+1 152 aubb(nzsbb)=au(j) 153 au(j)=0.d0 ! QUESTION why??? --> erasing from Kfe which is now Kii 154 irowbb(nzsbb)=lbb 155 lbb=lbb+1 156 cycle loop 157 endif 158 enddo 159 endif 160 enddo loop 161 adbb(kbb)=ad(i) ! copy diag at bb 162 ad(i)=1.d0 ! setting as 1 in Kfe, (to use at Kii) 163 kbb=kbb+1 164 enddo 165 jqbb(neqtot+1)=nzsbb+1 166 jqbi(neq(1)+1)=nzsbi+1 167 jqib(neqtot+1)=nzsib+1 ! must be ++1 jq: tot subdiags +1 168! 169 170 do i=1,neqtot 171 icolbb(i) = jqbb(i+1) - jqbb(i) 172 enddo 173 174 return 175 end 176