fe_matrix.c

矩阵运算〉〉 〉〉 〉〉 〉〉 。

 711:                   switch( gmatrix->eRep ) {
 712:                     case  INDIRECT :
 713:                         gmatrix->uMatrix.daa[k-1][m-1] += Ge->uMatrix.daa[j-1][l-1];
 714:                         break;
 715:                     case  SKYLINE :
 716:                         if( l >= j ) {
 717:                             iMin = (int) MIN( k, m );
 718:                             iMax = (int) MAX( k, m );
 719:                             if( (iMax-iMin+1) <= gmatrix->uMatrix.daa[iMax-1][0] )
 720:                                 gmatrix->uMatrix.daa[iMax-1][iMax-iMin+1]
 721:                                 += Ge->uMatrix.daa[j-1][l-1];
 722:                             else
 723:                                 gmatrix->uMatrix.daa[iMax-1][iMax-iMin+1]
 724:                                 =  Ge->uMatrix.daa[j-1][l-1];
 725:                         }
 726:                         break;
 727:                     default:
 728:                         FatalError("In Assemble_Global() : Undefined spA->eRep",
 729:                                     (char *) NULL);
 730:                         break;
 731:                   }
 732:                   
 733:                   if( UNITS_SWITCH == ON ) {
 734:                      UnitsCopy( &(gmatrix->spColUnits[m-1]), &(Ge->spColUnits[l-1]) );
 735:                      UnitsCopy( &(gmatrix->spRowUnits[k-1]), &(Ge->spRowUnits[j-1]) );
 736:                   }
 737:                }
 738:            }
 739: #ifdef DEBUG
 740:            printf(" local dof : j = %d \t global dof : k = %d \n", j,k); 
 741: #endif
 742:         }
 743:     }
 744: 
 745: #ifdef DEBUG
 746:        printf("*** Leave Assemble_Global()\n");
 747: #endif
 748: 
 749:     return(gmatrix);
 750: }
 751: 
 752: 
 753: 
 754: /* 
 755:  *  -----------------------------------------------
 756:  *  Assemble Global Load: Equiv Nodal Load  
 757:  *  -----------------------------------------------
 758:  */ 
 759: 
 760: #ifdef __STDC__
 761: MATRIX *Assemble_Global_Load(MATRIX *gmatrix, EFRAME *frame, int *Ld, MATRIX *Ge)
 762: #else
 763: MATRIX *Assemble_Global_Load(gmatrix, frame, Ld, Ge)
 764: MATRIX *gmatrix;
 765: EFRAME *frame;
 766: int    *Ld;
 767: MATRIX *Ge;
 768: #endif
 769: {
 770: ELEMENT          *el;
 771: int j,l,k,m, size_ld;
 772: int           length;
 773: int     UNITS_SWITCH;
 774: 
 775:     size_ld = Ld[0];
 776:     UNITS_SWITCH = CheckUnits();
 777: 
 778: #ifdef DEBUG
 779:        printf("*** Enter Assemble_Global_Load() : length of Ld = %3d : ", size_ld);
 780: #endif
 781: 
 782:     for(j = 1; j <= size_ld; j++) {
 783:         k = Ld[j];
 784:         if(k > 0) {
 785:           gmatrix->uMatrix.daa[k-1][0] += Ge->uMatrix.daa[j-1][0];
 786:           if( UNITS_SWITCH == ON ) {
 787:               UnitsCopy( &(gmatrix->spRowUnits[k-1]),  &(Ge->spRowUnits[j-1]) );
 788:           }
 789:         }
 790:     }
 791:     if( UNITS_SWITCH == ON )
 792:        ZeroUnits( &(gmatrix->spColUnits[0]) );
 793: 
 794: #ifdef DEBUG
 795:        printf("*** Leave Assemble_Global_Load()\n");
 796: #endif
 797: 
 798:     return(gmatrix);
 799: }
 800: 
 801: /* ======================================================= */
 802: /* Set Element Matrix - Stiffness/Mass for Feap problem    */
 803: /* Input  - frame pointer, element number,p array          */
 804: /* Output - Local Element Stiffness/Mass                   */
 805: /* ======================================================= */
 806: 
 807: #ifdef __STDC__
 808: MATRIX *Element_Matrix(ARRAY *p, int isw)
 809: #else
 810: MATRIX *Element_Matrix(p, isw)
 811: ARRAY    *p;
 812: int     isw;
 813: #endif
 814: {
 815:    p = elmlib(p,isw);
 816:    return(p->stiff);
 817: }
 818: 
 819: /* ========================================================== */
 820: /* Set Element Equivalent Load -                              */
 821: /* Input  - p working array , task id :isw = EQUIV_NODAL_LOAD */
 822: /* Output - Local Element Equivalent nodal load               */
 823: /* ========================================================== */
 824: 
 825: #ifdef __STDC__
 826: MATRIX *Element_Equiv(ARRAY *p, int isw)
 827: #else
 828: MATRIX *Element_Equiv(p, isw)
 829: ARRAY    *p;
 830: int     isw;
 831: #endif
 832: {
 833:    p = elmlib(p,isw);
 834:    return(p->equiv_nodal_load);
 835: }
 836: 
 837: /* ======================================================= */
 838: /* Set Element Vector -  Lumped Mass for Feap problem      */
 839: /* Input  - frame pointer, element number,p array          */
 840: /* Output - Local Lumped  Mass                             */
 841: /* ======================================================= */
 842: 
 843: #ifdef __STDC__
 844: QUANTITY *Element_Vector(ARRAY *p, int isw)
 845: #else
 846: QUANTITY *Element_Vector(p, isw)
 847: ARRAY     *p;
 848: int      isw;
 849: #endif
 850: {
 851: 
 852:   /*printf("*** In Element_VECTOR. Flag 1. Elem_# = %d\n",p->elmt_no); */
 853:   p = elmlib(p,isw);
 854:   return(p->nodal_loads);
 855: }
 856: 
 857: /* ===================================================================== */
 858: /* Rigidbody_conn                                                        */
 859: /* Input - frame pointer frp, array ptr checks if any node of an element */
 860: /*         is connected to a rigid body. It labels it by adding          */
 861: /*         a value of (one thousand * rigid_body no) to such nodes       */
 862: /*                                                                       */
 863: /* NOTE : Nodes attached to rigid body are scaled by a factor 1,000,000  */
 864: /*        in the working-element p->ix[] data structure.                 */
 865: /* ===================================================================== */
 866: 
 867: #ifdef __STDC__
 868: int Rigidbody_conn(EFRAME *frp, ARRAY *p)
 869: #else
 870: int Rigidbody_conn(frp,p)
 871: EFRAME *frp;
 872: ARRAY  *p;
 873: #endif
 874: {
 875: int i,j,k,nn, no_node,out;
 876: out = NO;
 877: 
 878:     for(i = 1; i<=frp->no_rigid;i++) {
 879:         no_node = frp->rigid[i-1].nodes_conn;
 880: 
 881:         for(j = 1; j<=no_node;j++){
 882:             nn = frp->rigid[i-1].in[j];
 883:             for(k = 1; k<=p->nodes_per_elmt;k++) {
 884:                 if(p->node_connect[k] == nn){
 885:                    p->node_connect[k] = p->node_connect[k] + i * 1000000;
 886:                    frp->node[nn -1].rb_num = i;
 887:                    out = YES;
 888:                 }
 889:             }
 890:         }
 891:     }
 892: 
 893:     return(out);
 894: }
 895: 
 896: 
 897: /* =================================================================== */
 898: /* Transform Stiff Matrix of member connected to Rigid Body            */
 899: /* transforms action at j of member to working pointt f of Rigid Body. */
 900: /*                                                                     */
 901: /* Transforms element S. Matrix Ke to T*Ke*Tt                          */
 902: /* Input - EFRAME *frp;  ARRAY *p;                                     */
 903: /*         Member Stiff Matrix MS                                      */
 904: /* Output- Modified Member Stiff Matrix MS                             */
 905: /* =================================================================== */
 906: 
 907: #ifdef __STDC__
 908: MATRIX *Transform_Stiff_Matrix(EFRAME *frp, ARRAY *p, MATRIX *MS)
 909: #else
 910: MATRIX *Transform_Stiff_Matrix(frp, p, MS)
 911: EFRAME     *frp;
 912: ARRAY        *p;
 913: MATRIX      *MS;
 914: #endif
 915: {
 916: int    nn,i,j,k,rb_no,node_no,rb_ty;
 917: double **Tfj, **Ti, **TiT, **ktit;
 918: double cgx,cgy,cgz;
 919: 
 920:     Tfj = (double **)MatrixAllocIndirectDouble(p->dof_per_node,p->dof_per_node);
 921:     Ti  = (double **)MatrixAllocIndirectDouble(p->size_of_stiff,p->size_of_stiff);
 922:     TiT = (double **)MatrixAllocIndirectDouble(p->size_of_stiff,p->size_of_stiff);
 923: 
 924:     for(i = 1; i<=p->nodes_per_elmt; i++) { 
 925:         nn = p->node_connect[i];
 926:         if(nn > 1000000) {
 927:            node_no = nn - rb_no * 1000000;
 928:            rb_no = nn/1000000; /* find rigid body number */
 929: 
 930:     /* ---- FIX LATER 
 931:            if(frp->ndm == 3) {
 932:                if(frp->basic_elmt_type == TRUSS_3D)
 933:                  rb_ty = TRUSS_3D;
 934:               else
 935:                  rb_ty = FRAME_3D;
 936:             }
 937:            else 
 938:                  rb_ty = frp->rigid[rb_no -1].rb_type; 
 939:            */ 
 940: 
 941:            cgx = p->coord[1][i].value - frp->rigid[rb_no -1].xcg.value;
 942:            cgy = p->coord[2][i].value - frp->rigid[rb_no -1].ycg.value;
 943: 
 944: /******** fix later
 945:            if(frp->ndm >= 3)
 946:               cgz = p->coord[3][i] - frp->rigid[rb_no -1].zcg;
 947:            else
 948:               cgz = 0.0;
 949: ****************************/
 950: 
 951:            Tfj = (double **) Transformation_Matrix(Tfj,rb_ty, p->dof_per_node,cgx,cgy,cgz);
 952:            Tfj = (double **)       Modify_T_Matrix(Tfj,  frp, p->dof_per_node, rb_no);
 953:        }
 954:        else 
 955:          Tfj = Transformation_Matrix(Tfj,rb_ty, p->dof_per_node,0.0,0.0,0.0);
 956: 
 957:        /* get other terms by symm          */
 958:        /* printf("node at end = %d\n", i); */
 959: 
 960:        for(j = 1; j<=p->dof_per_node; j++) {
 961:           for(k = 1; k<=p->dof_per_node; k++) {
 962:                Ti[(i-1)*p->dof_per_node+j][ (i-1)*p->dof_per_node+ k] = Tfj[j][k];
 963:               TiT[(i-1)*p->dof_per_node+k][ (i-1)*p->dof_per_node+ j] = Tfj[j][k];
 964:           }
 965:        }
 966:     }
 967: 
 968:    /* ------------------ */
 969:    /* Transformed Matrix */
 970:    /* ------------------ */
 971: 
 972:    ktit = (double **) dMatrixMult(MS->uMatrix.daa,p->size_of_stiff,p->size_of_stiff, TiT, p->size_of_stiff, p->size_of_stiff);
 973:    MS->uMatrix.daa   = (double **) dMatrixMultRep(MS->uMatrix.daa, Ti, p->size_of_stiff, p->size_of_stiff, ktit, p->size_of_stiff, p->size_of_stiff);
 974: 
 975:    MatrixFreeIndirectDouble(Tfj, p->dof_per_node);
 976:    MatrixFreeIndirectDouble(Ti, p->size_of_stiff);
 977:    MatrixFreeIndirectDouble(TiT, p->size_of_stiff);
 978:    MatrixFreeIndirectDouble(ktit, p->size_of_stiff);
 979: 
 980:    return(MS);
 981: }
 982: 
 983: 
 984: /* ============================================================ */
 985: /* Transform Mass Matrix of RigidBody:                          */
 986: /* Transforms R body Mass Matrix M =p->s to Mp = Tpc.M.Ttpc.    */
 987: /* transforms action from                                       */
 988: /*  mass centres   to  working point of Rigid Body              */
 989: /*                                                              */
 990: /*  here pcx,pcy,pcz are cordinates of mass centre c of RBODY   */
 991: /*                       wrt working point p of RBODY           */
 992: /* Input -       frame pointer frp, Array ptr p.                */
 993: /*                Member Stiff Matrix MS                        */
 994: /* Output-       Modified Member Stiff Matrix MS                */
 995: /* ============================================================ */
 996: 
 997: #ifdef __STDC__
 998: MATRIX *Transform_Rigid_Body_Mass_Matrix(EFRAME *frp, ARRAY *p, MATRIX *MS)
 999: #else
1000: MATRIX *Transform_Rigid_Body_Mass_Matrix(frp, p, MS)
1001: EFRAME   *frp;
1002: ARRAY    *p;
1003: MATRIX   *MS;
1004: #endif
1005: {
1006: int nn, i,j,k,rb_no,node_no,rb_ty;
1007: double **Tpc, **Ttpc,**MTtpc;
1008: double pcx,pcy,pcz;
1009: 
1010:      /*----------------------------------------------------------------------*/
1011:      /* Alloc_Matrix                                                         */
1012:      /*----------------------------------------------------------------------*/
1013: 
1014:         Tpc  = MatrixAllocIndirectDouble(p->dof_per_node,p->dof_per_node);
1015:         Ttpc = MatrixAllocIndirectDouble(p->dof_per_node,p->dof_per_node);
1016: 
1017:      /*----------------------------------------------------------------------*/
1018:      /* find rigid body number */
1019:      /*----------------------------------------------------------------------*/
1020: 
1021:      rb_no = p->elmt_no;         
1022:      rb_ty = frp->rigid[rb_no -1].rb_type;
1023: 
1024:      /*----------------------------------------------------------------------*/
1025:      /* Coordinates of mass centre c  wrt working point p(geom centre)       */
1026:      /*    = global coord of cg c - global coord of wkg pt  p                */ 
1027:      /*----------------------------------------------------------------------*/
1028: 
1029:      pcx =     p->work_section[4].value - frp->rigid[rb_no -1].xcg.value;
1030:      pcy =     p->work_section[5].value - frp->rigid[rb_no -1].ycg.value;
1031:      pcz = 0.0;
1032: /***********
1033:      if( frp->ndm >= 3)
1034:         pcz =    p->work_section[6] - frp->rigid[rb_no -1].zcg;
1035: ***************/
1036:      Tpc = (double **) Transformation_Matrix(Tpc,rb_ty, p->dof_per_node,pcx,pcy,pcz);
1037: 
1038:      /*----------------------------------------------------------------------*/
1039:      /* Modify depending on the dofs constrained   */
1040:      /*----------------------------------------------------------------------*/
1041: 
1042:      Tpc = (double **) Modify_T_Matrix(Tpc, frp,p->dof_per_node,rb_no);
1043: 
1044:     /* Transpose of Tpc matrix */
1045: 
1046:      for(j = 1; j<=p->dof_per_node;j++) 
1047:        for(k = 1; k<=p->dof_per_node;k++) 
1048:           Ttpc[j][k] = Tpc[k][j];
1049:       
1050:       /* Product : M.Ttpc           */
1051:       /* Product : Mp = Tpc. M.Ttpc */
1052: 
1053:       MTtpc = (double **) dMatrixMult(MS->uMatrix.daa, p->dof_per_node, p->dof_per_node, Ttpc, p->dof_per_node, p->dof_per_node);
1054:       MS->uMatrix.daa    = (double **) dMatrixMultRep(MS->uMatrix.daa, Tpc, p->dof_per_node, p->dof_per_node, MTtpc, p->dof_per_node, p->dof_per_node);
1055: 
1056:       MatrixFreeIndirectDouble(Tpc, p->dof_per_node);
1057:       MatrixFreeIndirectDouble(Ttpc, p->dof_per_node);
1058:       MatrixFreeIndirectDouble(MTtpc, p->dof_per_node);
1059: 
1060:       return(MS);
1061: }
1062: 
1063: 
1064: /* ======================================================*/
1065: /* Transformation matrix Tfj (=Tpj) for point j to f(=p) */