elmt_fiber2d.c

有限元程序

            if( isw == LOAD_MATRIX ) {
               Ltrans  = MatrixTranspose( L );
               temp_m1 = MatrixMult( Ltrans, Q );   /* element nodal force in global coord. */
               MatrixFree( Ltrans );
            }

            if( isw == STRESS ) {
               R = Rigid_Body_Rotation_2d( elmt_length );
               Rtrans = MatrixTranspose( R );
               temp_m1 = MatrixMult( Rtrans, Q );   /* element nodal force in local coord. */
               MatrixFree( R );
               MatrixFree( Rtrans );
            }

            MatrixFree( L );

            /* Assign force values */

            for( ii=0 ; ii < p->size_of_stiff ; ++ii )
                p->nodal_loads[ii].value = temp_m1->uMatrix.daa[ii][0];

            if( p->elmt_state == 0 ) {
                MatrixFree( K );
                MatrixFree( Q );
                MatrixFree( q );
            }

            MatrixFree( temp_m1 );

            /* Assign force units */

            if( UNITS_SWITCH == ON ) {
                SetUnitsOn();
                switch( UnitsType ) {
                   case SI:
                   case SI_US:
                        dp_force  = DefaultUnits("N");
                        dp_length = DefaultUnits("m");
                        break;
                   case US:
                        dp_force  = DefaultUnits("lbf");
                        dp_length = DefaultUnits("in");
                        break;
                }

                /* node 1 */

                UnitsCopy( p->nodal_loads[0].dimen, dp_force );
                UnitsCopy( p->nodal_loads[1].dimen, dp_force );
                UnitsMultRep( p->nodal_loads[2].dimen, dp_force, dp_length );

                /* node 2 */

                UnitsCopy( p->nodal_loads[3].dimen,
                           p->nodal_loads[0].dimen );
                UnitsCopy( p->nodal_loads[4].dimen,
                           p->nodal_loads[1].dimen );
                UnitsCopy( p->nodal_loads[5].dimen,
                           p->nodal_loads[2].dimen );

                free((char *) dp_force->units_name);
                free((char *) dp_length->units_name);
                free((char *) dp_force);
                free((char *) dp_length);
            }

            if(isw == STRESS ) {
               xx  = 0.5*(p->coord[0][0].value + p->coord[0][1].value);   /* xx = 0.5(x1+x2) */
               yy  = 0.5*(p->coord[1][0].value + p->coord[1][1].value);   /* yy = 0.5(y1+y2) */
               printf("\n");
               printf("Elmt No %3d : \n", p->elmt_no);
               switch( UNITS_SWITCH ) {
                  case ON :
                       printf("Coords (X,Y) = (%10.3f %s, %10.3f %s)\n", 
                          xx/elmt_length.dimen->scale_factor,elmt_length.dimen->units_name,
                          yy/elmt_length.dimen->scale_factor,elmt_length.dimen->units_name);
                       printf("\n");

                       /* node i */

                       printf(" Fx1 = %13.5e %s  Fy1 = %13.5e %s  Mz1 = %13.5e %s\n",
                          p->nodal_loads[0].value/p->nodal_loads[0].dimen->scale_factor,
                          p->nodal_loads[0].dimen->units_name,
                          p->nodal_loads[1].value/p->nodal_loads[1].dimen->scale_factor,
                          p->nodal_loads[1].dimen->units_name,
                          p->nodal_loads[2].value/p->nodal_loads[2].dimen->scale_factor,
                          p->nodal_loads[2].dimen->units_name);

                       /* node j */

                       printf(" Fx2 = %13.5e %s  Fy2 = %13.5e %s  Mz2 = %13.5e %s\n",
                          p->nodal_loads[3].value/p->nodal_loads[3].dimen->scale_factor,
                          p->nodal_loads[3].dimen->units_name,
                          p->nodal_loads[4].value/p->nodal_loads[4].dimen->scale_factor,
                          p->nodal_loads[4].dimen->units_name,
                          p->nodal_loads[5].value/p->nodal_loads[5].dimen->scale_factor,
                          p->nodal_loads[5].dimen->units_name);
                          printf("\n");

                      /* Member Forces */

                      printf(" Axial Force : x-direction = %13.5e %s\n",
                         -p->nodal_loads[0].value/p->nodal_loads[0].dimen->scale_factor,
                          p->nodal_loads[0].dimen->units_name);
                      printf(" Shear Force : y-direction = %13.5e %s\n",
                          p->nodal_loads[1].value/p->nodal_loads[1].dimen->scale_factor,
                          p->nodal_loads[1].dimen->units_name);
                      printf("\n");
                      break;
                 case OFF :
                      printf("Coords (X,Y) = (%10.3f , %10.3f )\n", xx, yy);
                      printf("\n");

                      /* node i */

                      printf(" Fx1 = %13.5e   Fy1 = %13.5e   Mz1 = %13.5e \n",
                         p->nodal_loads[0].value, 
                         p->nodal_loads[1].value, 
                         p->nodal_loads[2].value);

                      /* node j */

                      printf(" Fx2 = %13.5e   Fy2 = %13.5e   Mz2 = %13.5e \n",
                         p->nodal_loads[3].value,
                         p->nodal_loads[4].value,
                         p->nodal_loads[5].value);
                      printf("\n");

                      /* Member Forces */

                      printf(" Axial Force : x-direction = %13.5e \n", -p->nodal_loads[0].value);
                      printf(" Shear Force : y-direction = %13.5e \n",  p->nodal_loads[1].value);
                      printf("\n");
                      break;
                  default:
                      break;
               }
            }
            break;  /* end of case STRESS, LOAD_MATRIX */

       case STRESS_MATRIX:

            /* save element nodal forces in working array */

            if( UNITS_SWITCH == ON )
                SetUnitsOff();

            /* --------------------------------------------------------- */
            /* Two cases to deal with:                                   */
            /*                                                           */
            /* p->elmt_state == 0 : ElmtStateDet() has not been called.  */
            /*    Therefore, element stresses and loads must be computed */
            /*    from scratch.                                          */
            /* p->elmt_state == 1 : ElmtStateDet() has been called.      */
            /*    Therefore, element stresses and loads are saveed and   */
            /*    only need to be retrieved from Aladdin database.       */
            /* --------------------------------------------------------- */

            if( p->elmt_state == 1 ) {
                L = Element_Transformation_2d( p->coord, elmt_length );
                Q = p->Q_saved;
            }

            if( p->elmt_state == 0 ) { 

                fiber     = p->fiber_ptr->fiber;
                abscissas = (double *)MyCalloc( no_section, sizeof(double) );
                weights   = (double *)MyCalloc( no_section, sizeof(double) );
                Gauss_Lobatto( abscissas, weights, no_integ_pt );

                /* Allocate memory for element force/stiffness */

                kx = MatrixAllocIndirect( "kx", DOUBLE_ARRAY, 3, 3 );
                bx = MatrixAllocIndirect( "bx", DOUBLE_ARRAY, 3, 3 );

                /* Construct element flexibility matrix by looping over sections */

                for( sec=0 ; sec < no_section; ++sec ) {
                   xi    = elmt_length.value/2.0*(abscissas[sec]+1);
                   scale = weights[sec]*elmt_length.value/2.0;

                   Force_Interpolation_Matrix_2d( bx, xi, elmt_length );
                   Section_Tangent_Stiffness_2d( kx, fiber, no_fiber, no_shear, sec, elmt_no );

                   fx = MatrixInverse( kx );

                   bxtrans = MatrixTranspose( bx );
                   temp_m1 = MatrixMult( bxtrans, fx );
                   temp_m2 = MatrixMult( temp_m1, bx );

                   MatrixFree( fx );
                   MatrixFree( temp_m1 );
                   MatrixFree( bxtrans );

                   if( sec == 0 )
                      F = MatrixScale( temp_m2, scale );
                   else {
                      temp_m1 = MatrixScale( temp_m2, scale );
                      MatrixAddReplace( F, temp_m1 );
                      MatrixFree( temp_m1 );
                   }

                   MatrixFree( temp_m2 );
                }

                /* Invert element flexibility to get element-level stiffness */

                K = MatrixInverse( F );

                /* Clean up ..... */

                free((char *) abscissas);
                free((char *) weights);
                MatrixFree( kx );
                MatrixFree( bx );
                MatrixFree( F );

                /* --------------------------------------------------------- */
                /* Put displacement matrix into one column format i.e.,      */
                /* Transfer                                                  */
                /*                                                           */
                /*  p->displ=[ px1, px2;   to temp_m1 = pe = [ px1;          */
                /*             py1, py2;                       py1;          */
                /*             pz1, pz2 ]                      pz1;          */
                /*                                             px2;          */
                /*                                             py2;          */
                /*                                             pz2 ]         */
                /* --------------------------------------------------------- */

                temp_m1  = MatrixAllocIndirect( (char *)NULL, DOUBLE_ARRAY, 6, 1 );
                for( ii=0 ; ii < p->dof_per_node ; ++ii )
                for( jj=0 ; jj < p->nodes_per_elmt ; ++jj )
                   temp_m1->uMatrix.daa[ii+jj*p->dof_per_node][0] = p->displ->uMatrix.daa[ii][jj];

                L = Element_Transformation_2d( p->coord, elmt_length );
                q = MatrixMult( L, temp_m1 );
                Q = MatrixMult( K, q );
                MatrixFree( temp_m1 );
            }

            /* Compute nodal forces in local element coordinates */

            R       = Rigid_Body_Rotation_2d( elmt_length );
            Rtrans  = MatrixTranspose( R );
            temp_m1 = MatrixMult( Rtrans, Q );

            /* Assign force values */

            for( ii=0 ; ii < p->size_of_stiff ; ++ii )
                 p->nodal_loads[ii].value = temp_m1->uMatrix.daa[ii][0];

            /* Cleanup memory from working arrays */

            if( p->elmt_state == 0 ) {
                MatrixFree( K );
                MatrixFree( Q );
                MatrixFree( q );
            }

            MatrixFree( R );
            MatrixFree( Rtrans );
            MatrixFree( L );
            MatrixFree( temp_m1 );

            /* Assign force units */

            if( UNITS_SWITCH == ON ) {
               SetUnitsOn();
               if( CheckUnitsType() == SI ) {
                   dp_length = DefaultUnits("m");
                   dp_force  = DefaultUnits("N");
               } else if (CheckUnitsType() == US ) {
                   dp_length = DefaultUnits("in");
                   dp_force  = DefaultUnits("lbf");
               } 

               /* node 1 */

               UnitsCopy( p->nodal_loads[0].dimen, dp_force );
               UnitsCopy( p->nodal_loads[1].dimen, dp_force );
               UnitsMultRep( p->nodal_loads[2].dimen, dp_force, dp_length );

               /* node 2 */

               UnitsCopy( p->nodal_loads[3].dimen, p->nodal_loads[0].dimen );
               UnitsCopy( p->nodal_loads[4].dimen, p->nodal_loads[1].dimen );
               UnitsCopy( p->nodal_loads[5].dimen, p->nodal_loads[2].dimen );

               /* Transfer nodal coordinates and forces/moment to working array */

               UnitsCopy( &(p->stress->spColUnits[0]), dp_length );
               UnitsCopy( &(p->stress->spColUnits[1]), dp_length );

               UnitsCopy( &(p->stress->spColUnits[2]), dp_force  );
               UnitsCopy( &(p->stress->spColUnits[3]), dp_force  );
               UnitsMultRep( &(p->stress->spColUnits[4]), dp_force, dp_length );

               /* Zero out units buffer */

               ZeroUnits( &(p->stress->spRowUnits[0]) );
               ZeroUnits( &(p->stress->spRowUnits[1]) );

               /* Release working memory */

               free((char *) dp_force->units_name);
               free((char *) dp_length->units_name);
               free((char *) dp_force);
               free((char *) dp_length);
            }

            /* Transfer coordinates to working stress matrix */

            p->stress->uMatrix.daa[0][0] = p->coord[0][0].value;
            p->stress->uMatrix.daa[0][1] = p->coord[1][0].value;
            p->stress->uMatrix.daa[1][0] = p->coord[0][1].value;
            p->stress->uMatrix.daa[1][1] = p->coord[1][1].value;

            /* Transfer internal loads to working stress matrix */

            p->stress->uMatrix.daa[0][2] = p->nodal_loads[0].value;
            p->stress->uMatrix.daa[0][3] = p->nodal_loads[1].value;
            p->stress->uMatrix.daa[0][4] = p->nodal_loads[2].value;
            p->stress->uMatrix.daa[1][2] = p->nodal_loads[3].value;
            p->stress->uMatrix.daa[1][3] = p->nodal_loads[4].value;
            p->stress->uMatrix.daa[1][4] = p->nodal_loads[5].value;


#ifdef DEBUG
       printf("*** Finished case STRESS_MATRIX : isw = %d\n", isw );

       for (ii = 1; ii <= 2; ii = ii + 1 ) {
            for (jj = 1; jj <= 5; jj = jj + 1 ) {
                 printf(" %10.4f ", p->stress->uMatrix.daa[ ii-1 ][ jj-1 ]);
            }
            printf("\n");
       }

       /* MatrixPrintVar ( p->stress ); */
       printf("*** Leaving STRESS_MATRIX case\n" );
#endif

            break; 
       case MASS_MATRIX:   /* compute element mass matrix */

          /* mbar = density * area */

          if( p->work_section[6].value != 0.0 )    /* mbar = weight/gravity */
             mbar = p->work_section[6].value/9.80665;
          else if( (p->work_material[5].value > 0.0) && (p->work_section[10].value > 0.0) )
             mbar = p->work_material[5].value * p->work_section[10].value;
          else
             FatalError("\nError in input: Need density value to calculate mass matrix\n",(char *)NULL);

          cs = p->coord[0][1].value - p->coord[0][0].value;
          sn = p->coord[1][1].value - p->coord[1][0].value;
          cs = cs / elmt_length.value;
          sn = sn / elmt_length.value;

          switch( p->type ) {  /* mass type */
             case LUMPED :
                /* use lumped mass matrix of FRAME_2D element */
                p->stiff = beamms( p, p->stiff, p->type, mbar, elmt_length.value, cs, sn,
                                   p->size_of_stiff, p->dof_per_node );
                break;

             case CONSISTENT :
             default :