📄 fe_matrix.c
字号:
356: *
357: * Usage : load = Form_Internal_Load(displ)
358: * or load = Form_Internal_Load(displ, displ_incr)
359: *
360: * Boundary Info : INPUT -- bound_id = 0, no restraint,
361: * bound_id > 0, restrained, displ specified
362: * bound_id is changed based on INPUT bound_id
363: * : OUTPUT -- bound_id < 0, restrained displ specified
364: * bound_id > 0, no restraint
365: *
366: * ==============================================================
367: */
368:
369: #ifdef __STDC__
370: MATRIX *Form_Internal_Load( MATRIX *displ, ... )
371: #else
372: MATRIX *Form_Internal_Load(va_alist)
373: va_dcl
374: #endif
375: {
376: #ifndef __STDC__
377: MATRIX *displ;
378: #endif
379:
380: va_list arg_ptr;
381: MATRIX *displ_incr;
382: MATRIX *load;
383: ELEMENT *ep;
384: int node_no, elmt_no;
385: int i, j, k, jj, *Ld;
386: int elmt_attr_no;
387: int length1, length2;
388: int UNITS_SWITCH;
389: int itemp;
390:
391: #ifdef DEBUG
392: printf("*** Enter Form_Internal_Load()\n");
393: #endif
394:
395: UNITS_SWITCH = CheckUnits();
396:
397: /* [a] : read input */
398:
399: #ifdef __STDC__
400: va_start(arg_ptr, displ);
401: #else
402: va_start(arg_ptr);
403: displ = va_arg(arg_ptr, MATRIX *);
404: #endif
405:
406: if(array->material_name != NULL && strcmp(array->material_name, "ELASTIC")) {
407: displ_incr = va_arg(arg_ptr, MATRIX *);
408: } else
409: displ_incr = (MATRIX *) NULL;
410:
411: va_end(arg_ptr);
412:
413: /* [b] : Allocate memory for internal load vector */
414:
415: load = MatrixAllocIndirect("Internal Load", DOUBLE_ARRAY, TNEQ, 1);
416:
417: /* [c] : Assemble internal load vector from element level loads */
418:
419: for(elmt_no = 1; elmt_no <= frame->no_elements; elmt_no++) {
420:
421: ep = &frame->element[elmt_no-1];
422: if(displ_incr == (MATRIX *) NULL)
423: ep->esp->state = 0; /* Calculate internal load as elastic behaviour */
424: else
425: ep->esp->state = 1; /* Need check for ELASTIC PLASTIC STATE */
426:
427: array = Assign_p_Array(frame, elmt_no, array, STRESS);
428: array = elmlib(array, PROPTY);
429:
430: /* Transfer fixed displacements */
431:
432: for(i=1; i <= array->nodes_per_elmt; i++) {
433:
434: k = 1;
435: node_no = ep->node_connect[i-1];
436: for(j = 1; j <= frame->no_dof; j++) {
437:
438: switch((int) array->nodes_per_elmt) {
439: case 2: case 3:
440: jj = frame->node[node_no - 1].bound_id[j-1];
441: if(jj > 0) {
442:
443: /* displacement */
444:
445: array->displ->uMatrix.daa[j-1][i-1] = displ->uMatrix.daa[jj-1][0];
446: if( UNITS_SWITCH == ON ) {
447: UnitsCopy(&(array->displ->spRowUnits[j-1]), &(displ->spRowUnits[jj-1]));
448: UnitsCopy(&(array->displ->spColUnits[i-1]), &(displ->spColUnits[0]));
449: }
450:
451: /* displacement increments */
452:
453: if(displ_incr != (MATRIX *) NULL) {
454: array->displ_incr->uMatrix.daa[j-1][i-1] = displ_incr->uMatrix.daa[jj-1][0];
455: if(UNITS_SWITCH == ON) {
456: UnitsCopy(&(array->displ_incr->spRowUnits[j-1]), &(displ->spRowUnits[jj-1]));
457: UnitsCopy(&(array->displ_incr->spColUnits[i-1]), &(displ->spColUnits[0]));
458: }
459: }
460: } else {
461:
462: array->displ->uMatrix.daa[j-1][i-1] = frame->node[node_no -1].disp[j-1].value;
463:
464: if(UNITS_SWITCH == ON ) {
465:
466: /* displacement */
467:
468: UnitsCopy(&(array->displ->spRowUnits[j-1]), frame->node[node_no -1].disp[j-1].dimen);
469: UnitsCopy(&(array->displ->spColUnits[i-1]), &(displ->spColUnits[0]));
470:
471: /* displcement increments */
472:
473: UnitsCopy(&(array->displ_incr->spRowUnits[j-1]), frame->node[node_no -1].disp[j-1].dimen);
474: UnitsCopy(&(array->displ_incr->spColUnits[i-1]), &(displ->spColUnits[0]));
475: }
476: }
477: break;
478: case 4: case 8:
479: jj = frame->node[node_no - 1].bound_id[j-1];
480: if(jj > 0) {
481:
482: /* displacement */
483:
484: array->displ->uMatrix.daa[k-1][i-1] = displ->uMatrix.daa[jj-1][0];
485:
486: if( UNITS_SWITCH == ON ) {
487: UnitsCopy(&(array->displ->spRowUnits[k-1]), &(displ->spRowUnits[jj-1]));
488: UnitsCopy(&(array->displ->spColUnits[i-1]), &(displ->spColUnits[0]));
489: }
490:
491: /* incremental displacement */
492:
493: if(displ_incr != (MATRIX *) NULL) {
494:
495: array->displ_incr->uMatrix.daa[k-1][i-1] = displ_incr->uMatrix.daa[jj-1][0];
496:
497: if(UNITS_SWITCH == ON) {
498: UnitsCopy(&(array->displ_incr->spRowUnits[k-1]), &(displ->spRowUnits[jj-1]));
499: UnitsCopy(&(array->displ_incr->spColUnits[i-1]), &(displ->spColUnits[0]));
500: }
501: }
502: } else {
503: array->displ->uMatrix.daa[k-1][i-1] = frame->node[node_no -1].disp[j-1].value;
504:
505: if(UNITS_SWITCH == ON ) {
506:
507: /* displacement */
508:
509: UnitsCopy(&(array->displ->spRowUnits[k-1]), frame->node[node_no -1].disp[j-1].dimen);
510: UnitsCopy(&(array->displ->spColUnits[i-1]), &(displ->spColUnits[0]));
511:
512: /* displcement increments */
513:
514: UnitsCopy(&(array->displ_incr->spRowUnits[k-1]), frame->node[node_no-1].disp[j-1].dimen);
515: UnitsCopy(&(array->displ_incr->spColUnits[i-1]), &(displ->spColUnits[0]));
516: }
517: }
518: k = k + 1;
519: break;
520: default:
521: break;
522: }
523: }
524: }
525:
526: array = elmlib(array, LOAD_MATRIX);
527: Ld = Destination_Array(frame, elmt_no);
528:
529:
530:
531: if( UNITS_SWITCH == ON )
532: ZeroUnits(&(load->spColUnits[0]));
533:
534: for(j = 1; j <= Ld[0]; j++) {
535: k = Ld[j];
536: if(k > 0) {
537: load->uMatrix.daa[k-1][0] += array->nodal_loads[j-1].value;
538: if( UNITS_SWITCH == ON )
539: UnitsCopy(&(load->spRowUnits[k-1]), array->nodal_loads[j-1].dimen);
540: }
541: }
542:
543: #ifdef DEBUG
544: printf("DIAGNOSTIC : elmt node %4d\n", elmt_no );
545: dMatrixPrint("elmt diplacements", array->displ->uMatrix.daa, 3,2);
546: printf("destination :");
547: for(itemp = 1; itemp <= Ld[0]; itemp++) {
548: printf(" %10d : ", Ld[itemp]);
549: }
550: printf("\n");
551: printf("elmt nodal loads :");
552: for(itemp = 1; itemp <= Ld[0]; itemp++) {
553: printf(" %10.2e : ", array->nodal_loads[itemp-1].value);
554: }
555: printf("\n");
556: dMatrixPrint("Internal Load", load->uMatrix.daa, TNEQ, 1);
557: #endif
558:
559: free((char *) Ld);
560: }
561:
562: #ifdef DEBUG
563: printf("*** Leave Form_Internal_Load()\n");
564: #endif
565:
566: return(load);
567: }
568:
569:
570: /*
571: * =================================================
572: * UTILITY Functions for Stiffness and Mass matrices
573: * =================================================
574: */
575:
576: /*
577: * -----------------------------------------
578: * [a] : Transfer Destination Array to *Ld[]
579: * -----------------------------------------
580: */
581:
582: #ifdef __STDC__
583: int *Destination_Array(EFRAME *frame, int elmt_no)
584: #else
585: int *Destination_Array(frame, elmt_no)
586: EFRAME *frame;
587: int elmt_no;
588: #endif
589: {
590: ELEMENT *el;
591: int i,no_dof_per_elmt, *Ld;
592:
593: #ifdef DEBUG
594: printf(" enter Destination_Array() \n");
595: printf(" in Destination_Array() : \n");
596: printf(" : frame->no_dof = %d \n", frame->no_dof);
597: printf(" : frame->no_nodes_per_elmt = %d \n", frame->no_nodes_per_elmt);
598: #endif
599:
600: no_dof_per_elmt = frame->no_dof*frame->no_nodes_per_elmt;
601:
602: #ifdef DEBUG
603: printf(" : no_dof_per_elmt = %d \n", no_dof_per_elmt);
604: #endif
605:
606: Ld = iVectorAlloc((no_dof_per_elmt + 1));
607:
608: el = &frame->element[elmt_no - 1];
609: for(i = 1; i <= no_dof_per_elmt + 1; i++) {
610: Ld[i-1] = el->d_array[i-1];
611: }
612:
613: #ifdef DEBUG
614: printf(" leaving Destination_Array() \n");
615: #endif
616: return(Ld);
617: }
618:
619: /*
620: * ----------------------------------------------------
621: * [c] : Compute Destination Array *Ld[] for Rigid Body
622: * ----------------------------------------------------
623: */
624:
625: #ifdef __STDC__
626: int *Destination_Array_for_Rigid_Body(EFRAME *frp, ARRAY *p, int elmt_no, int dflg)
627: #else
628: int *Destination_Array_for_Rigid_Body(frp, p, elmt_no, dflg)
629: EFRAME *frp;
630: ARRAY *p;
631: int elmt_no, dflg;
632: #endif
633: {
634: ELEMENT *el;
635: int no_nodes_per_elmt, no_dof;
636: int i,j,k,jj,ii,iid,mate_no, rb_ty;
637: int *vld;
638:
639: no_nodes_per_elmt = p->nodes_per_elmt; /* MIN(p->nodes_per_elmt, frp->no_nodes_per_elmt);*/
640: no_dof = p->dof_per_node; /* MIN(p->dof_per_node,frp->no_dof);*/
641: vld = iVectorAlloc(no_dof*no_nodes_per_elmt);
642:
643: rb_ty = frp->rigid[elmt_no -1].rb_type ;
644: ii = frp->rigid[elmt_no -1].in[1] ; /* node no of node conn by rbody for dofs */
645:
646: for(j=1;j<=no_dof;j++){ /* loop over dofs */
647: jj = frp->rigid[elmt_no -1].rest_dof[j];
648: if(jj >= 1 && rb_ty != 7) /* dofs not connected */
649: vld[j] = 0; /*k;*/
650: else { /* jj=0 i:e; rest dofs & for rb = 7; storey type rbody */
651: k = frp->node[ii - 1].bound_id[j]; /* equation nos assigned to dofs */
652: vld[j] = k;
653: }
654: }
655:
656: vld[0] = no_dof*no_nodes_per_elmt;
657:
658: return(vld);
659: }
660:
661:
662: /*
663: * ------------------------------------------------------------------
664: * Assemble Global Matrix : Stiffness and Mass
665: * gmatrix : Global matrix TNEQxTNEQ
666: * Ge : Local matrix array->size_of_stiffxarray->size_of_stiff
667: * ------------------------------------------------------------------
668: */
669:
670: #ifdef __STDC__
671: MATRIX *Assemble_Global(MATRIX *gmatrix, EFRAME *frame, int *Ld, MATRIX *Ge)
672: #else
673: MATRIX *Assemble_Global(gmatrix, frame, Ld, Ge)
674: MATRIX *gmatrix;
675: EFRAME *frame;
676: int *Ld;
677: MATRIX *Ge;
678: #endif
679: {
680: ELEMENT *el;
681: int j,l,k,m, size_ld;
682: int length1, length2;
683: int iMin, iMax;
684: int UNITS_SWITCH;
685: int linked, count, ik;
686:
687: size_ld = Ld[0];
688: UNITS_SWITCH = CheckUnits();
689:
690: #ifdef DEBUG
691: printf("*** Enter Assemble_Global() : length of Ld = %3d : ", size_ld);
692: for(j = 1; j <= size_ld; j++) {
693: printf(" %3d ", Ld[j]);
694: }
695: printf("\n");
696: #endif
697:
698: for(j = 1; j <= size_ld; j++) {
699: k = Ld[j];
700:
701: count = 0;
702: for(ik = 1; ik <= size_ld; ik++) {
703: if(k == Ld[ik]) count = count + 1;
704: }
705:
706: linked = (count > 1);
707: if(k > 0 && linked == 0) {
708: for(l = 1; l <= size_ld; l++){
709: m = Ld[l];
710: if(m > 0){
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -