fe_matrix.c

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

1421: ARRAY    *p;
1422: QUANTITY *fvector, *pl;
1423: int      elmt_no;
1424: #endif
1425: {
1426: int      j,k;
1427: int      i,iid,ii;
1428: ELEMENT  *el;
1429: int      *ldofg;
1430: char     *el_type;
1431: int      elmt_attr_no;
1432: 
1433:     el           = &frp->element[elmt_no -1];
1434:     elmt_attr_no = el->elmt_attr_no;
1435:     el_type      = frp->eattr[elmt_attr_no-1].elmt_type;
1436:     ldofg        = frp->eattr[elmt_no-1].map_ldof_to_gdof;
1437: 
1438:     for(i=1;i<=p->nodes_per_elmt;i++){
1439:         ii  = el->node_connect[i];
1440:         iid = (ii-1) * frp->no_dof;
1441:         for(j = 1; j<= p->dof_per_node; j++){
1442:             k = ldofg[j-1];
1443:             if(k > 0)
1444:               fvector[iid + k].value += pl[(i-1)*p->dof_per_node+ j].value;
1445:         }
1446:     }
1447: 
1448:     return(fvector);
1449: }
1450: 
1451: /* =============================== */
1452: /* Function to Assemble Nodal Load */
1453: /* =============================== */
1454: 
1455: #ifdef __STDC__
1456: QUANTITY *Assemble_Nodal_Load(EFRAME *frp, QUANTITY *fv)
1457: #else
1458: QUANTITY *Assemble_Nodal_Load(frp, fv)
1459: EFRAME   *frp;
1460: QUANTITY *fv;
1461: #endif
1462: {
1463: NODE               *np;
1464: NODE_LOADS        *nlp;
1465: int i,node, dof, nload;
1466: 
1467:  /* Initialize load vector */
1468: 
1469:     for(dof=1; dof<=frp->no_eq; dof++)
1470:         fv[dof].value = 0.;
1471: 
1472:     for(nload = 1; nload <= frp->no_node_loads; nload++) {
1473:         nlp  = &frp->nforces[nload-1];
1474:         node = nlp->node_f;
1475:         for(i =1; i <= frp->no_dof; i++)
1476:             fv[(node -1) *  frp->no_dof + i ].value = nlp->fn[i].value;
1477:     }
1478: 
1479:     return (fv);
1480: }
1481: 
1482: /* ======================================= */
1483: /*   Function to Assemble  Gravity  Load   */
1484: /* ======================================= */
1485: 
1486: #ifdef __STDC__
1487: QUANTITY *Assemble_Gravity_Load(EFRAME *frp, QUANTITY *fv)
1488: #else
1489: QUANTITY *Assemble_Gravity_Load(frp, fv)
1490: EFRAME   *frp;
1491: QUANTITY *fv;
1492: #endif
1493: {
1494: NODE               *np;
1495: NODE_LOADS        *nlp;
1496: int i,node, dof, nload;
1497: 
1498:    /* FILL IN DETAILS LATER */
1499: 
1500:    return (fv);
1501: }
1502: 
1503: /* =========================================== */
1504: /*   Function to Apply the Boundary conditions */
1505: /*   to the global stiffness matrix            */
1506: /* =========================================== */
1507: 
1508: #ifdef __STDC__
1509: double **Boundary_Conditions(EFRAME *frp)
1510: #else
1511: double **Boundary_Conditions(frp)
1512: EFRAME *frp;
1513: #endif
1514: {
1515: NODE                *np;
1516: NODE_LOADS         *nlp;
1517: int  i,node, dof, nload;
1518: 
1519:    for(node=1;node<=frp->no_nodes;node++) {
1520:        np   = &frp->node[node-1];
1521:        for(i =1; i<= frp->no_dof; i++) {
1522:        if(np->bound_id[i] >0) {
1523:            for(dof=1;dof<=TDOF;dof++) {
1524:                K[dof][(node-1)*frp->no_dof +i] = 0.;
1525:                K[(node-1)*frp->no_dof + 1][dof] = 0.;
1526:            }
1527: 
1528:            K[(node-1)*frp->no_dof + 1][(node-1)*frp->no_dof+ 1] = 0.;
1529:            F[(node-1)*frp->no_dof + 1]   = 0.;
1530:        }
1531:        }
1532:    }
1533: 
1534:    return (K);
1535: }
1536: 
1537: /* ==================================== */
1538: /*   Function to compute nodal forces   */
1539: /* ==================================== */
1540: 
1541: #ifdef __STDC__
1542: QUANTITY *Nodal_Forces(EFRAME *frp)
1543: #else
1544: QUANTITY *Nodal_Forces(frp)
1545: EFRAME   *frp;
1546: #endif
1547: {
1548: int i,j;
1549: 
1550:     /* Initialize nodal force vector : Compute [K]*{d} = {Fn} */
1551: 
1552:      for(i=1; i<=TDOF; i++) Fn[i].value = 0.;
1553:          for(i=1;i<=TDOF;i++)
1554:          for(j=1;j<=TDOF;j++)
1555:              Fn[i].value += K[i][j]*d[j];
1556: 
1557:      return (Fn);
1558: }
1559: 
1560: /* =======================================  */
1561: /*   Function to Assemble  Centrifugal  Load*/
1562: /* =======================================  */
1563: 
1564: #ifdef __STDC__
1565: QUANTITY *Assemble_Ctrfgl_Load(EFRAME *frp, QUANTITY *fv)
1566: #else
1567: QUANTITY *Assemble_Ctrfgl_Load(frp,fv)
1568: EFRAME   *frp;
1569: QUANTITY *fv;
1570: #endif
1571: {
1572: NODE               *np;
1573: NODE_LOADS        *nlp;
1574: int i,node, dof, nload;
1575: 
1576:  /* to be changed */
1577: 
1578:    for(dof=1; dof<=frp->no_eq; dof++)
1579:       fv[dof].value   = 0.;
1580: 
1581:    for(nload = 1; nload <= frp->no_node_loads; nload++) {
1582: 
1583:      nlp  = &frp->nforces[nload-1];
1584:      node = nlp->node_f;
1585:          for(i =1; i <= frp->no_dof; i++)
1586:              fv[(node -1) *  frp->no_dof + i ].value = nlp->fn[i].value;
1587: 
1588:    }
1589: 
1590:    return (fv);
1591: }
1592: 
1593: 
1594: 
1595: /* Functions for design rule checking */
1596: 
1597: #ifdef __STDC__
1598: MATRIX *Get_Coord(MATRIX *m)
1599: #else
1600: MATRIX *Get_Coord(m)
1601: MATRIX *m;
1602: #endif
1603: {
1604: MATRIX *coord;
1605: int    nodeno;
1606: int        ii;
1607: 
1608: #ifdef DEBUG
1609:        printf("*** Enter Get_Coord()\n");
1610: #endif
1611: 
1612:        nodeno = (int) m->uMatrix.daa[0][0];
1613:        coord  = MatrixAllocIndirect("Node Coord", DOUBLE_ARRAY, 1, frame->no_dimen);
1614: 
1615:        if( CheckUnits() == ON ) {
1616:           for( ii=0 ; ii<frame->no_dimen ; ii++ ) {
1617:              coord->uMatrix.daa[0][ii] = frame->node[nodeno-1].coord[ii].value;
1618:              UnitsCopy( &(coord->spColUnits[ii]), frame->node[nodeno-1].coord[ii].dimen );
1619:           }
1620:           ZeroUnits( &(coord->spRowUnits[0]) );
1621:        }
1622:        else {
1623:           for( ii=0 ; ii<frame->no_dimen ; ii++ )
1624:              coord->uMatrix.daa[0][ii] = frame->node[nodeno-1].coord[ii].value;
1625:        }
1626: 
1627: #ifdef DEBUG
1628:        printf("*** Leave Get_Coord()\n");
1629: #endif
1630: 
1631:     return(coord);
1632: }
1633: 
1634: #ifdef __STDC__
1635: MATRIX *Get_Node(MATRIX *m)
1636: #else
1637: MATRIX *Get_Node(m)
1638: MATRIX *m;
1639: #endif
1640: {
1641: MATRIX *connect;
1642: int      elmtno;
1643: int          ii;
1644: 
1645: #ifdef DEBUG
1646:        printf("*** Enter Get_Node()\n");
1647: #endif
1648: 
1649:        elmtno   = (int) m->uMatrix.daa[0][0];
1650:        connect  = MatrixAllocIndirect("Node Connect", DOUBLE_ARRAY, 1, frame->no_nodes_per_elmt);
1651: 
1652:        if( CheckUnits() == ON ) {
1653:           for( ii=0 ; ii<frame->no_nodes_per_elmt ; ii++ ) {
1654:              connect->uMatrix.daa[0][ii] = (double) frame->element[elmtno-1].node_connect[ii];
1655:              ZeroUnits( &(connect->spColUnits[ii]) );
1656:           }
1657:           ZeroUnits( &(connect->spRowUnits[0]) );
1658:        }
1659:        else {
1660:           for( ii=0 ; ii<frame->no_nodes_per_elmt ; ii++ )
1661:              connect->uMatrix.daa[0][ii] = (double) frame->element[elmtno-1].node_connect[ii];
1662:        }
1663: 
1664: #ifdef DEBUG
1665:        printf("*** Leave Get_Node()\n");
1666: #endif
1667: 
1668:     return(connect);
1669: }
1670: 
1671: #ifdef __STDC__
1672: MATRIX *Get_Displ(MATRIX *m1, MATRIX *m2)
1673: #else
1674: MATRIX *Get_Displ(m1,m2)
1675: MATRIX *m1, *m2;
1676: #endif
1677: {
1678: MATRIX *displ;
1679: int    nodeno;
1680: int     ii,jj;
1681: int UNITS_SWITCH;
1682: 
1683: #ifdef DEBUG
1684:        printf("*** Enter Get_Displ()\n");
1685: #endif
1686: 
1687:        nodeno = (int) m1->uMatrix.daa[0][0];
1688:        displ  = MatrixAllocIndirect("Node Displ", DOUBLE_ARRAY, 1, frame->no_dof);
1689:        UNITS_SWITCH = CheckUnits();
1690: 
1691:           switch( UNITS_SWITCH ) {
1692:             case ON:
1693:                    for( ii=0 ; ii<frame->no_dof ; ii++ ) {
1694:                       jj = frame->node[nodeno-1].bound_id[ii];
1695:                       if(jj > 0) {
1696:                          if(m2->spRowUnits[jj-1].units_name != NULL ) {
1697:                             UnitsTypeConvert(&(m2->spRowUnits[jj-1]), UNITS_TYPE);
1698:                          }
1699:                          RadUnitsSimplify( &(m2->spRowUnits[jj-1]) );
1700:                          frame->node[nodeno-1].disp[ii].value = m2->uMatrix.daa[jj-1][0];
1701:                          UnitsCopy( frame->node[nodeno-1].disp[ii].dimen, &(m2->spRowUnits[jj-1]) );
1702:                       }
1703:                       else
1704:                          RadUnitsSimplify( frame->node[nodeno-1].disp[ii].dimen );
1705:                   }
1706:               break;
1707:            case OFF:
1708:                   for( ii=0 ; ii<frame->no_dof ; ii++){
1709:                      jj = frame->node[nodeno-1].bound_id[ii];
1710:                      if(jj > 0)
1711:                         frame->node[nodeno-1].disp[ii].value = m2->uMatrix.daa[jj-1][0];
1712:                   }
1713:               break;
1714:            default:
1715:               break;
1716:           }
1717: 
1718:        if( UNITS_SWITCH == ON ) {
1719:           for( ii=0 ; ii<frame->no_dof ; ii++ ) {
1720:              displ->uMatrix.daa[0][ii] = frame->node[nodeno-1].disp[ii].value;
1721:              UnitsCopy( &(displ->spColUnits[ii]), frame->node[nodeno-1].disp[ii].dimen );
1722:           }
1723:           ZeroUnits( &(displ->spRowUnits[0]) );
1724:        }
1725:        else {
1726:           for( ii=0 ; ii<frame->no_dof ; ii++ )
1727:              displ->uMatrix.daa[0][ii] = frame->node[nodeno-1].disp[ii].value;
1728:        }
1729: 
1730: #ifdef DEBUG
1731:        printf("*** Leave Get_Displ()\n");
1732: #endif
1733: 
1734:     return(displ);
1735: }
1736: 
1737: #ifdef __STDC__
1738: MATRIX *Get_Stress(MATRIX *m1, MATRIX *m2)
1739: #else
1740: MATRIX *Get_Stress(m1,m2)
1741: MATRIX *m1, *m2;
1742: #endif
1743: {
1744: MATRIX *stress;
1745: int    elmt_no;
1746: int UNITS_SWITCH;
1747: 
1748: ELEMENT        *ep;
1749: ELEMENT_ATTR   *eap;
1750: int node_no, elmt_attr_no;
1751: int  i,j,k,ii,jj;
1752: 
1753: #ifdef DEBUG
1754:        printf("*** Enter Get_Stress()\n");
1755: #endif
1756: 
1757:        elmt_no = (int) m1->uMatrix.daa[0][0];
1758:        stress  = MatrixAllocIndirect("Element Stress", DOUBLE_ARRAY, frame->no_nodes_per_elmt, frame->no_dof);
1759:        UNITS_SWITCH = CheckUnits();
1760: 
1761: #ifdef DEBUG
1762:         printf(" elmt No = %d \n", elmt_no);
1763: #endif
1764: 
1765:            array = Assign_p_Array(frame, elmt_no, array, STRESS);
1766:            array = elmlib(array, PROPTY);
1767: 
1768:            /* Transfer Fixed Displacements */
1769: 
1770:            ep           = &frame->element[elmt_no-1];
1771:            elmt_attr_no = ep->elmt_attr_no;  
1772: