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 subroutine mafillplhs(kon,ipkon,lakon,ne, 20 & ipompc,nodempc,coefmpc,nmpc,nactdoh,icolp,jqp, 21 & irowp,neqp,nzlp,nzsp,adbp,aubp, 22 & ipvar,var) 23! 24! filling the lhs pressure matrix in sparse matrix format 25! 26! it is assumed that the temperature MPC's also apply to the 27! pressure. Temperature MPC's are not allowed to contain 28! other variables than the temperature. 29! 30 implicit none 31! 32 character*8 lakon(*) 33! 34 integer kon(*),ipompc(*),nodempc(3,*),icolp(*),jqp(*),nzsp, 35 & nactdoh(*),konl(20),irowp(*),ipkon(*),ipvar(*),ne,nmpc,neqp, 36 & nzlp,i,j,jj,ll,id,id1,id2,ist,ist1,ist2,index,jdof1,jdof2, 37 & idof1,idof2,mpc1,mpc2,index1,index2,node1,node2,indexe,nope, 38 & i0 39! 40 real*8 coefmpc(*),sm(8,8),adbp(*),aubp(*),var(*),value 41! 42 i0=0 43! 44! determining nzlp 45! 46 nzlp=0 47 do i=neqp,1,-1 48 if(icolp(i).gt.0) then 49 nzlp=i 50 exit 51 endif 52 enddo 53! 54 do i=1,neqp 55 adbp(i)=0.d0 56 enddo 57 do i=1,nzsp 58 aubp(i)=0.d0 59 enddo 60! 61! loop over all fluid elements 62! 63 do i=1,ne 64! 65 if(ipkon(i).lt.0) cycle 66 if(lakon(i)(1:1).ne.'F') cycle 67 indexe=ipkon(i) 68 if(lakon(i)(4:4).eq.'8') then 69 nope=8 70 elseif(lakon(i)(4:4).eq.'4') then 71 nope=4 72 elseif(lakon(i)(4:4).eq.'6') then 73 nope=6 74 endif 75! 76 do j=1,nope 77 konl(j)=kon(indexe+j) 78 enddo 79! 80 call e_c3d_plhs(lakon(i),sm,i,ipvar,var) 81! 82 do jj=1,nope 83! 84 node1=kon(indexe+jj) 85 jdof1=nactdoh(node1) 86! 87 do ll=jj,nope 88! 89 node2=kon(indexe+ll) 90 jdof2=nactdoh(node2) 91! 92! check whether one of the DOF belongs to a SPC or MPC 93! 94 if((jdof1.gt.0).and.(jdof2.gt.0)) then 95 call add_sm_fl(aubp,adbp,jqp,irowp,jdof1,jdof2, 96 & sm(jj,ll),jj,ll) 97 elseif((jdof1.gt.0).or.(jdof2.gt.0)) then 98! 99! idof1: genuine DOF 100! idof2: nominal DOF of the SPC/MPC 101! 102 if(jdof1.le.0) then 103 idof1=jdof2 104 idof2=jdof1 105 else 106 idof1=jdof1 107 idof2=jdof2 108 endif 109 if(nmpc.gt.0) then 110 if(idof2.ne.2*(idof2/2)) then 111! 112! regular DOF / MPC 113! 114 id=(-idof2+1)/2 115 ist=ipompc(id) 116 index=nodempc(3,ist) 117 if(index.eq.0) cycle 118 do 119 idof2=nactdoh(nodempc(1,index)) 120 if(idof2.gt.0) then 121 value=-coefmpc(index)*sm(jj,ll)/ 122 & coefmpc(ist) 123 if(idof1.eq.idof2) value=2.d0*value 124 call add_sm_fl(aubp,adbp,jqp,irowp, 125 & idof1,idof2,value,i0,i0) 126 endif 127 index=nodempc(3,index) 128 if(index.eq.0) exit 129 enddo 130 cycle 131 endif 132 endif 133 else 134 idof1=jdof1 135 idof2=jdof2 136 mpc1=0 137 mpc2=0 138 if(nmpc.gt.0) then 139 if(idof1.ne.2*(idof1/2)) mpc1=1 140 if(idof2.ne.2*(idof2/2)) mpc2=1 141 endif 142 if((mpc1.eq.1).and.(mpc2.eq.1)) then 143 id1=(-idof1+1)/2 144 id2=(-idof2+1)/2 145 if(id1.eq.id2) then 146! 147! MPC id1 / MPC id1 148! 149 ist=ipompc(id1) 150 index1=nodempc(3,ist) 151 if(index1.eq.0) cycle 152 do 153 idof1=nactdoh(nodempc(1,index1)) 154 index2=index1 155 do 156 idof2=nactdoh(nodempc(1,index2)) 157 if((idof1.gt.0).and.(idof2.gt.0)) then 158 value=coefmpc(index1)*coefmpc(index2)* 159 & sm(jj,ll)/coefmpc(ist)/coefmpc(ist) 160 call add_sm_fl(aubp,adbp,jqp, 161 & irowp,idof1,idof2,value,i0,i0) 162 endif 163! 164 index2=nodempc(3,index2) 165 if(index2.eq.0) exit 166 enddo 167 index1=nodempc(3,index1) 168 if(index1.eq.0) exit 169 enddo 170 else 171! 172! MPC id1 / MPC id2 173! 174 ist1=ipompc(id1) 175 index1=nodempc(3,ist1) 176 if(index1.eq.0) cycle 177 do 178 idof1=nactdoh(nodempc(1,index1)) 179 ist2=ipompc(id2) 180 index2=nodempc(3,ist2) 181 if(index2.eq.0) then 182 index1=nodempc(3,index1) 183 if(index1.eq.0) then 184 exit 185 else 186 cycle 187 endif 188 endif 189 do 190 idof2=nactdoh(nodempc(1,index2)) 191 if((idof1.gt.0).and.(idof2.gt.0)) then 192 value=coefmpc(index1)*coefmpc(index2)* 193 & sm(jj,ll)/coefmpc(ist1)/coefmpc(ist2) 194 if(idof1.eq.idof2) value=2.d0*value 195 call add_sm_fl(aubp,adbp,jqp, 196 & irowp,idof1,idof2,value,i0,i0) 197 endif 198! 199 index2=nodempc(3,index2) 200 if(index2.eq.0) exit 201 enddo 202 index1=nodempc(3,index1) 203 if(index1.eq.0) exit 204 enddo 205 endif 206 endif 207 endif 208 enddo 209 enddo 210 enddo 211! 212 return 213 end 214