root/src/dfdbj.c

/* [<][>][^][v][top][bottom][index][help] */

DEFINITIONS

This source file includes following definitions.
  1. dfdbj

   1 /*
   2      CalculiX - A 3-dimensional finite element program
   3               Copyright (C) 1998-2015 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 /*                                                                       */
  20 /*     author: Reinhold Fischer                                          */
  21 
  22 #include <stdio.h>
  23 #include <math.h>
  24 #include <stdlib.h>
  25 #include <string.h>
  26 #include "CalculiX.h"
  27 
  28 void dfdbj(double *bcont,double **dbcontp,ITG *neq,ITG *nope,ITG *konl,
  29            ITG* nactdof,double *s,double *z,ITG *ikmpc,ITG *ilmpc,
  30            ITG *ipompc,ITG *nodempc,ITG *nmpc,double *coefmpc,
  31            double *fnl,ITG *nev,ITG **ikactcontp,ITG **ilactcontp,
  32            ITG *nactcont,ITG *nactcont_,ITG *mi, ITG *cyclicsymmetry,
  33            ITG *izdof, ITG *nzdof){
  34 
  35   ITG j,j1,jdof,kdof,k,k1,l,id,index,ist,id1,ist1,index1,id2,ist2,index2,
  36       jdbcontcol,i1,i3,i4,mt=mi[1]+1,im,*ikactcont=*ikactcontp,
  37       *ilactcont=*ilactcontp,kdofm1;
  38 
  39   double d1,sl,*dbcont=*dbcontp;
  40   
  41   for(j=0; j<*nope; j++){
  42       i1=mt*(konl[j]-1)+1;
  43       for(j1=0; j1<3; j1++){
  44           jdof=nactdof[i1+j1];
  45           if(jdof!=0){
  46               jdof--;
  47               FORTRAN(nident,(ikactcont,&jdof,nactcont,&id));
  48               do{
  49                   if(id>0){
  50                       if(ikactcont[id-1]==jdof){
  51                           jdbcontcol=ilactcont[id-1];
  52                           break;
  53                       }
  54                   }
  55                   (*nactcont)++;
  56                   if(*nactcont>*nactcont_){
  57                       *nactcont_=(ITG)(1.1**nactcont_);
  58                       RENEW(ikactcont,ITG,*nactcont_);
  59                       RENEW(ilactcont,ITG,*nactcont_);
  60                       RENEW(dbcont,double,*nev**nactcont_);
  61                   }
  62                   k=*nactcont-1;
  63                   l=k-1;
  64                   while(k>id){
  65                       ikactcont[k]=ikactcont[l];
  66                       ilactcont[k--]=ilactcont[l--];
  67                   }
  68                   jdbcontcol=*nactcont;
  69                   ikactcont[id]=jdof;
  70                   ilactcont[id]=*nactcont;
  71 //                memset(&dbcont[(*nactcont-1)**nev],0,sizeof(double)**nev);
  72                   DMEMSET(dbcont,(*nactcont-1)**nev,*nactcont**nev,0.);
  73                   break;
  74               }while(1);
  75               bcont[jdof]-=fnl[j*3+j1];
  76               i4=(jdbcontcol-1)**nev;
  77               i3=(3*j+j1);
  78               for(k=0; k<*nope; k++){
  79                   for(k1=0; k1<3; k1++){
  80                       sl=s[(3*k+k1)*100+i3];
  81                       kdof=nactdof[mt*(konl[k]-1)+k1+1];
  82                       if(kdof!=0){
  83                           if(!(*cyclicsymmetry)){
  84                           for(l=0; l<*nev; l++){
  85                               dbcont[i4+l]-=sl*z[(long long)l**neq+kdof-1];
  86                           }
  87                         }else{
  88                           kdofm1=kdof-1;
  89                           FORTRAN(nident,(izdof,&kdofm1,nzdof,&id));
  90                           if(id!=0){
  91                             if(izdof[id-1]==kdofm1){
  92                               for(l=0; l<*nev; l++){
  93                                 dbcont[i4+l]-=sl*z[l**nzdof+id-1];
  94                               }
  95                             }else{printf("*ERROR in dfdbj\n");FORTRAN(stop,());}
  96                           }else{printf("*ERROR in dfdbj\n");FORTRAN(stop,());}
  97                         }
  98                       }
  99                       else{
 100                           kdof=8*(konl[k]-1)+k1+1;
 101                           FORTRAN(nident,(ikmpc,&kdof,nmpc,&id));
 102                           if(id>0){
 103                               id--;
 104                               if(ikmpc[id]==kdof){
 105                                   id=ilmpc[id];
 106                                   ist=ipompc[id-1];
 107                                   ist--;
 108                                   index=nodempc[ist*3+2];
 109                                   if(index==0) continue;
 110                                   index--;
 111                                   do{
 112                                       kdof=nactdof[mt*(nodempc[index*3]-1)+nodempc[index*3+1]];
 113                                       d1=sl*coefmpc[index]/coefmpc[ist];
 114                                       if(kdof!=0){
 115                                           if(!(*cyclicsymmetry)){
 116                                           for(l=0; l<*nev; l++){
 117                                             dbcont[i4+l]+=d1*z[(long long)l**neq+kdof-1];
 118                                           }
 119                                         }
 120                                       }else{
 121                                         kdofm1=kdof-1;
 122                                         FORTRAN(nident,(izdof,&kdofm1,nzdof,&id));
 123                                         if(id!=0){
 124                                           if(izdof[id-1]==kdofm1){
 125                                             for(l=0; l<*nev; l++){
 126                                               dbcont[i4+l]+=d1*z[l**nzdof+id-1];
 127                                             }
 128                                           }else{printf("*ERROR in dfdbj\n");FORTRAN(stop,());}
 129                                         }else{printf("*ERROR in dfdbj\n");FORTRAN(stop,());}
 130                                       }
 131                                       index=nodempc[index*3+2];
 132                                       if(index==0) break;
 133                                       index--;
 134                                   }while(1);
 135                               }
 136                           }
 137                       }
 138                   }
 139               }
 140           }
 141           else{
 142               jdof=8*(konl[j]-1)+j1+1;
 143               FORTRAN(nident,(ikmpc,&jdof,nmpc,&id1));
 144               if(id1>0){
 145                   id1--;
 146                   if(ikmpc[id1]==jdof){
 147                       id1=ilmpc[id1];
 148                       ist1=ipompc[id1-1];
 149                       ist1--;
 150                       index1=nodempc[ist1*3+2];
 151                       if(index1==0) continue;
 152                       index1--;
 153                       do{
 154                           jdof=nactdof[mt*(nodempc[index1*3]-1)+nodempc[index1*3+1]];
 155                           if(jdof!=0){
 156                               jdof--;
 157                               FORTRAN(nident,(ikactcont,&jdof,nactcont,&id));
 158                               do{
 159                                   if(id>0){
 160                                       if(ikactcont[id-1]==jdof){
 161                                           jdbcontcol=ilactcont[id-1];
 162                                       }
 163                                   }
 164                                   (*nactcont)++;
 165                                   if(*nactcont>*nactcont_){
 166                                       *nactcont_=(ITG)(1.1**nactcont_);
 167                                       RENEW(ikactcont,ITG,*nactcont_);
 168                                       RENEW(ilactcont,ITG,*nactcont_);
 169                                       RENEW(dbcont,double,*nev**nactcont_);
 170                                   }
 171                                   k=*nactcont-1;
 172                                   l=k-1;
 173                                   do{
 174                                       ikactcont[k]=ikactcont[l];
 175                                       ilactcont[k--]=ilactcont[l--];
 176                                   }while(k>id);
 177                                   jdbcontcol=*nactcont;
 178                                   ikactcont[id]=jdof;
 179                                   ilactcont[id]=*nactcont;
 180 //                                memset(&dbcont[(*nactcont-1)**nev],0,sizeof(double)**nev);              
 181                                   DMEMSET(dbcont,(*nactcont-1)**nev,*nactcont**nev,0.);
 182                                   break;
 183                               }while(1);
 184                               bcont[jdof]+=coefmpc[index1]*fnl[j*3+j1]/coefmpc[ist1];
 185                               i4=(jdbcontcol-1)**nev;
 186                               i3=(3*j+j1);
 187                               for(k=0; k<*nope; k++){
 188                                   for(k1=0; k1<3; k1++){
 189                                       sl=s[(3*k+k1)*100+i3];
 190                                       kdof=nactdof[mt*(konl[k]-1)+k1+1];
 191                                       if(kdof!=0){
 192                                           d1=sl*coefmpc[index1]/coefmpc[ist1];
 193                                           if(!(*cyclicsymmetry)){
 194                                             for(l=0; l<*nev; l++){
 195                                               dbcont[i4+l]+=d1*z[(long long)l**neq+kdof-1];
 196                                             }
 197                                           }else{
 198                                             kdofm1=kdof-1;
 199                                             FORTRAN(nident,(izdof,&kdofm1,nzdof,&id));
 200                                             if(id!=0){
 201                                               if(izdof[id-1]==kdofm1){
 202                                                 for(l=0; l<*nev; l++){
 203                                                   dbcont[i4+l]+=d1*z[l**nzdof+id-1];
 204                                                 }
 205                                               }else{printf("*ERROR in dfdbj\n");FORTRAN(stop,());}
 206                                             }else{printf("*ERROR in dfdbj\n");FORTRAN(stop,());}
 207                                           }
 208                                       }
 209                                       else{
 210                                           kdof=8*(konl[k]-1)+k1+1;
 211                                           FORTRAN(nident,(ikmpc,&kdof,nmpc,&id2));
 212                                           if(id2>0){
 213                                               id2--;
 214                                               if(ikmpc[id2]==kdof){
 215                                                   id2=ilmpc[id2];
 216                                                   ist2=ipompc[id2-1];
 217                                                   ist2--;
 218                                                   index2=nodempc[ist2*3+2];
 219                                                   if(index2==0) continue;
 220                                                   index2--;
 221                                                   do{
 222                                                       kdof=nactdof[mt*(nodempc[index2*3]-1)+nodempc[index2*3+1]];
 223                                                       if(kdof!=0){
 224                                                           d1=sl*coefmpc[index1]*coefmpc[index2]/(coefmpc[ist1]*coefmpc[ist2]);
 225                                                           if(!(*cyclicsymmetry)){
 226                                                             for(l=0; l<*nev; l++){
 227                                                               dbcont[i4+l]-=d1*z[(long long)l**neq+kdof-1];
 228                                                             }
 229                                                           }else{
 230                                                             kdofm1=kdof-1;
 231                                                             FORTRAN(nident,(izdof,&kdofm1,nzdof,&id));
 232                                                             if(id!=0){
 233                                                               if(izdof[id-1]==kdofm1){
 234                                                                 for(l=0; l<*nev; l++){
 235                                                                   dbcont[i4+l]-=d1*z[l**nzdof+id-1];
 236                                                                 }
 237                                                               }else{printf("*ERROR in dfdbj\n");FORTRAN(stop,());}
 238                                                             }else{printf("*ERROR in dfdbj\n");FORTRAN(stop,());}
 239                                                           }
 240                                                       }
 241                                                       index2=nodempc[index2*3+2];
 242                                                       if(index2==0) break;
 243                                                       index2--;
 244                                                   }while(1);
 245                                               }
 246                                           }
 247                                       }
 248                                   }
 249                               }
 250                           }
 251                           index1=nodempc[index1*3+2];
 252                           if(index1==0) break;
 253                           index1--;
 254                       }while(1);
 255                   }
 256               }
 257           }
 258       }
 259   }
 260   *dbcontp=dbcont;
 261   *ikactcontp=ikactcont;
 262   *ilactcontp=ilactcont;
 263 }
 264 
 265 /*!
 266   !     CalculiX - A 3-dimensional finite element program
 267   !              Copyright (C) 1998-2015 Guido Dhondt
 268   !
 269 !     This program is free software; you can redistribute it and/or
 270 !     modify it under the terms of the GNU General Public License as
 271 !     published by the Free Software Foundation(version 2);
 272 !     
 273 !
 274 !     This program is distributed in the hope that it will be useful,
 275 !     but WITHOUT ANY WARRANTY; without even the implied warranty of 
 276 !     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
 277 !     GNU General Public License for more details.
 278 !
 279 !     You should have received a copy of the GNU General Public License
 280 !     along with this program; if not, write to the Free Software
 281 !     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 282 !
 283         subroutine dfdbj(bcont,dbcont,neq,nope,konl,nactdof,s,z,
 284      &  ikmpc,ilmpc,ipompc,nodempc,nmpc,coefmpc,fnl,nev,iactcont,
 285      &  nactcont)
 286 !
 287 !     calculates the derivative of the contact forces with respect
 288 !     to the modal variables
 289 !
 290       implicit none
 291 !
 292       integer j,j1,neq,nope,konl(*),nactdof(0:3,*),jdof,kdof,
 293      &  k,k1,l,id,ikmpc(*),ilmpc(*),ipompc(*),nodempc(3,*),nmpc,
 294      &  index,ist,id1,ist1,index1,id2,ist2,index2,nev ,iactcont(*),
 295      &  nactcont,jdofcont
 296 !
 297       real*8 bcont(*),dbcont(nev,*),s(60,60),z(neq,*),coefmpc(*),
 298      &  fnl(3,9)
 299 !
 300       do j=1,nope
 301          do j1=1,3
 302             jdof=nactdof(j1,konl(j))
 303             if(jdof.ne.0) then
 304                call nident(iactcont,jdof,nactcont,id)
 305                jdofcont=0
 306                if(id.gt.0)then
 307                   if(iactcont(id).eq.jdof) then
 308                      jdofcont=id
 309                   endif
 310                endif
 311                if(jdofcont.eq.0) then
 312                   nactcont=nactcont+1
 313                   do k=nactcont,id+2,-1
 314                      iactcont(k)=iactcont(k-1)
 315                      do l=1,nev
 316                         dbcont(l,k)=dbcont(l,k-1)
 317                      enddo
 318                   enddo
 319                   jdofcont=id+1
 320                   iactcont(jdofcont)=jdof
 321                   do l=1,nev
 322                      dbcont(l,jdofcont)=0.d0
 323                   enddo
 324                endif
 325                bcont(jdof)=bcont(jdof)-fnl(j1,j)
 326                do k=1,nope
 327                   do k1=1,3
 328                      kdof=nactdof(k1,konl(k))
 329                      if(kdof.ne.0) then
 330                         do l=1,nev
 331                            dbcont(l,jdofcont)=dbcont(l,jdofcont)-
 332      &                       s(3*(j-1)+j1,3*(k-1)+k1)*z(kdof,l)
 333                         enddo
 334                      else
 335                         kdof=8*(konl(k)-1)+k1
 336                         call nident(ikmpc,kdof,nmpc,id)
 337                         if(id.gt.0) then
 338                            if(ikmpc(id).eq.kdof) then
 339                               id=ilmpc(id)
 340                               ist=ipompc(id)
 341                               index=nodempc(3,ist)
 342                               if(index.eq.0) cycle
 343                               do
 344                                  kdof=nactdof(nodempc(2,index),
 345      &                                        nodempc(1,index))
 346                                  if(kdof.ne.0) then
 347                                     do l=1,nev
 348                                        dbcont(l,jdofcont)=
 349      &                                      dbcont(l,jdofcont)+
 350      &                                   s(3*(j-1)+j1,3*(k-1)+k1)*
 351      &                                   coefmpc(index)*z(kdof,l)/
 352      &                                   coefmpc(ist)
 353                                     enddo
 354                                  endif
 355                                  index=nodempc(3,index)
 356                                  if(index.eq.0) exit
 357                               enddo
 358                            endif
 359                         endif
 360                      endif
 361                   enddo
 362                enddo
 363             else
 364                jdof=8*(konl(j)-1)+j1
 365                call nident(ikmpc,jdof,nmpc,id1)
 366                if(id1.gt.0) then
 367                   if(ikmpc(id1).eq.jdof) then
 368                      id1=ilmpc(id1)
 369                      ist1=ipompc(id1)
 370                      index1=nodempc(3,ist1)
 371                      if(index1.eq.0) cycle
 372                      do
 373                         jdof=nactdof(nodempc(2,index1),
 374      &                               nodempc(1,index1))
 375                         if(jdof.ne.0) then
 376                            call nident(iactcont,jdof,nactcont,id)
 377                            jdofcont=0
 378                            if(id.gt.0)then
 379                               if(iactcont(id).eq.jdof) then
 380                                  jdofcont=id
 381                               endif
 382                            endif
 383                            if(jdofcont.eq.0) then
 384                               nactcont=nactcont+1
 385                               do k=nactcont,id+2,-1
 386                                  iactcont(k)=iactcont(k-1)
 387                                  do l=1,nev
 388                                     dbcont(l,k)=dbcont(l,k-1)
 389                                  enddo
 390                               enddo
 391                               jdofcont=id+1
 392                               iactcont(jdofcont)=jdof
 393                               do l=1,nev
 394                                  dbcont(l,jdofcont)=0.d0
 395                               enddo
 396                            endif
 397                            bcont(jdofcont)=bcont(jdofcont)+
 398      &                          coefmpc(index1)*
 399      &                          fnl(j1,j)/coefmpc(ist1)
 400                            do k=1,nope
 401                               do k1=1,3
 402                                  kdof=nactdof(k1,konl(k))
 403                                  if(kdof.ne.0) then
 404                                     do l=1,nev
 405                                        dbcont(l,jdofcont)=
 406      &                                      dbcont(l,jdofcont)
 407      &                                   +s(3*(j-1)+j1,3*(k-1)+k1)
 408      &                                   *coefmpc(index1)*z(kdof,l)/
 409      &                                   coefmpc(ist1)
 410                                     enddo
 411                                  else
 412                                     kdof=8*(konl(k)-1)+k1
 413                                     call nident(ikmpc,kdof,nmpc,id2)
 414                                     if(id2.gt.0) then
 415                                        if(ikmpc(id2).eq.kdof) then
 416                                           id2=ilmpc(id2)
 417                                           ist2=ipompc(id2)
 418                                           index2=nodempc(3,ist2)
 419                                           if(index2.eq.0) cycle
 420                                           do
 421 !
 422 !                   translated to the left to avoid exceedance
 423 !                   of 72 columns
 424 !
 425                              kdof=nactdof(nodempc(2,index2),
 426      &                            nodempc(1,index2))
 427                              if(kdof.ne.0) then
 428                                 do l=1,nev
 429                                    dbcont(l,jdofcont)=dbcont(l,jdofcont)
 430      &                                  -s(3*(j-1)+j1,3*(k-1)+k1)
 431      &                                  *coefmpc(index1)
 432      &                                  *coefmpc(index2)*z(kdof,l)/
 433      &                                  (coefmpc(ist1)*coefmpc(ist2))
 434                                 enddo
 435                              endif
 436                              index2=nodempc(3,index2)
 437                              if(index2.eq.0) exit
 438 !
 439 !                   end of translation
 440 !
 441                                           enddo
 442                                        endif
 443                                     endif
 444                                  endif
 445                               enddo
 446                            enddo
 447                         endif
 448                         index1=nodempc(3,index1)
 449                         if(index1.eq.0) exit
 450                      enddo
 451                   endif
 452                endif
 453             endif
 454          enddo
 455       enddo
 456 !
 457       return
 458       end
 459   */

/* [<][>][^][v][top][bottom][index][help] */