FEATool  v1.8
Finite Element Analysis Toolbox
 All Files Functions Pages
lcoord_cell.m File Reference

Description

LCOORD_CELL Local cell coordinates for the degrees of freedom.

[ XI ] = LCOORD_CELL( NLDOF, N_SDIM, N_VC ) Computes the local cell coordinates for all local degree of freedom. NLDOF is an array specifying the number of local degrees of freedom on vertices, edges, faces, and cell interiors. N_SDIM is the number of space dimensions, and N_VC the number of vertices per cell.

Code listing

% Local cell vertex coordinates.
 if( n_vc==(n_sdim+1) )
   xiv = eye(n_sdim+1);
 else
   if( n_sdim==2 )
     xiv = [ -1  1 1 -1 ;
             -1 -1 1  1 ];
   elseif( n_sdim==3 )
     xiv = [ -1  1  1 -1 -1  1 1 -1 ;
             -1 -1  1  1 -1 -1 1  1 ;
             -1 -1 -1 -1  1  1 1  1 ];
   end
 end


 xi = [];

 for j=1:size(nLDof,2)
   dim = j-1;

   for i=1:size(nLDof,1)
     n_ldof_ij = nLDof(i,j);
     if( n_ldof_ij==0 )
       continue   % Skip.
     end

     switch( dim )

       case 0   % Vertices.

         xi = [ xi xiv ];

       case 1   % Edges.

         n_ec = n_vc*( n_sdim==2 ) + (6*(n_vc==4) + 12*(n_vc==8))*( n_sdim==3 );

         xie = [];
         for k_ldof=1:n_ldof_ij
           k_edge = mod(k_ldof-1,n_ec) + 1;

           dt  = ceil(k_ldof/n_ec) * 1/((n_ldof_ij/n_ec)+1);
           xie = [ xie lcoorde( k_edge, dt, n_vc, n_sdim ) ];
         end

         xi = [ xi xie ];

       case 2   % Faces.

         n_f_c = 4*(n_vc==4) + 6*(n_vc==8);

         if( n_vc==4 )

           switch( n_ldof_ij )
             case 4
               dt = ones(3,1)/3;
             case 12
               dt = [ 1/2 1/4 1/4 ;
                      1/4 1/2 1/4 ;
                      1/4 1/4 1/2 ];
             case 24
               dt = [ 3/5 1/5 1/5 ;
                      2/5 2/5 1/5 ;
                      1/5 3/5 1/5 ;
                      1/5 2/5 2/5 ;
                      1/5 1/5 3/5 ;
                      2/5 1/5 2/5 ]';
             otherwise
               warning(['lcoord_cell: dof positions for ',int2str(n_ldof_ij),' face dofs not supported on tetrahedra.'])
           end

         elseif( n_vc==8 )

           if( (n_ldof_ij/6)==4 )

             dt = [ 1 2 2 1; 1 1 2 2 ]/3;

           elseif( (n_ldof_ij/6)==9 )

             dt = [ 1 2 3 3 3 2 1 1 2; 1 1 1 2 3 3 3 2 2]/4;

           elseif( (n_ldof_ij/6)==16 )

             dt = [ 1 2 3 4 4 4 4 3 2 1 1 1 2 3 3 2; 1 1 1 1 2 3 4 4 4 4 3 2 2 2 3 3]/5;

           elseif( rem(sqrt(n_ldof_ij/6),1)==0 )

             xi1 = linspace( 0, 1, sqrt(n_ldof_ij/6) + 2 );
             [xi2,xi1] = meshgrid( xi1(2:end-1), xi1(2:end-1) );
             dt = [ xi1(:)'; xi2(:)' ];

           else
             warning(['lcoord_cell: dof positions for ',int2str(n_ldof_ij),' face dofs not supported on hexahedra.'])
           end

         end

         xif = [];
         for i_dt=1:size(dt,2)
           for k_f=1:n_f_c
             xif = [ xif lcoorde( k_f, dt(:,i_dt), n_vc, n_sdim ) ];
           end
         end

         xi = [ xi xif ];

       case 3   % Interior.

         xii = [];
         switch( n_sdim )

           case 1

             xi1 = linspace( 1, 0, n_ldof_ij+2 );
             xii = [ xi1(2:end-1); 1-xi1(2:end-1) ];

           case 2

             if( n_vc==3 )

               switch( n_ldof_ij )
                 case 1
                   xii = ones(3,1)/3;
                 case 3
                   xii = [ 1/2 1/4 1/4 ;
                           1/4 1/2 1/4 ;
                           1/4 1/4 1/2 ];
                 case 6
                   xii = [ 3/5 1/5 1/5 2/5 1/5 2/5 ;
                           1/5 3/5 1/5 2/5 2/5 1/5 ;
                           1/5 1/5 3/5 1/5 2/5 2/5 ];
                 otherwise
                   warning(['lcoord_cell: dof positions for ',int2str(n_ldof_ij),' interior face dofs not supported on triangles.'])
               end

             elseif( n_vc==4 )

               if( n_ldof_ij==3 )
                 xii = zeros(2,3);
               elseif( n_ldof_ij==4 )   % Use counterclockwise ordering for 4 internal quad dofs.
                 xii = [ -1 1 1 -1; -1 -1 1 1 ]/3;
               elseif( n_ldof_ij==9 )
                 xii = [-1  1 1 -1  0 1 0 -1 0;-1 -1 1  1 -1 0 1  0 0]/2;
               elseif( n_ldof_ij==16 )
                 xii = [-3/5  3/5 3/5 -3/5  -1/5  3/5 1/5 -3/5   1/5 3/5 -1/5 -3/5  -1/5  1/5 1/5 -1/5;
                        -3/5 -3/5 3/5  3/5  -3/5 -1/5 3/5  1/5  -3/5 1/5  3/5 -1/5  -1/5 -1/5 1/5  1/5];
               elseif( rem(sqrt(n_ldof_ij),1)==0 )
                 xi1 = linspace( -1, 1, sqrt(n_ldof_ij)+2 );
                 [xi2,xi1] = meshgrid( xi1(2:end-1), xi1(2:end-1) );
                 xii = [ xi1(:)'; xi2(:)' ];
               else
                 error(['lcoord_cell: quadrilateral cells with ',num2str(n_ldof_ij),' internal dofs not supported.'])
               end

             end

           case 3

             if( n_vc==4 )

               switch( n_ldof_ij )
                 case 1
                   xii = ones(4,1)/4;
                 case 4
                   xii = [ 2/5 1/5 1/5 1/5 ;
                           1/5 2/5 1/5 1/5 ;
                           1/5 1/5 2/5 1/5 ;
                           1/5 1/5 1/5 2/5 ];
                 otherwise
                   warning(['lcoord_cell: dof positions for ',int2str(n_ldof_ij),' interior face dofs not supported on tetrahedra.'])
               end

             elseif( n_vc==8 )

               if( n_ldof_ij==4 )
                 xii = zeros(3,4);
               elseif( n_ldof_ij==8 )   % Use counterclockwise ordering for 8 internal hex dofs.
                 xii = [ -1 1 1 -1 -1 1 1 -1; -1 -1 1 1 -1 -1 1 1; -1 -1 -1 -1 1 1 1 1 ]/3;
               elseif( n_ldof_ij==64 )
                 xii = [-3/5 -1/5  1/5  3/5 -3/5 -1/5  1/5  3/5 -3/5 -1/5  1/5  3/5 -3/5 -1/5  1/5  3/5 -3/5 -1/5  1/5  3/5 -3/5 -1/5  1/5  3/5 -3/5 -1/5  1/5  3/5 -3/5 -1/5  1/5  3/5 -3/5 -1/5  1/5  3/5 -3/5 -1/5  1/5  3/5 -3/5 -1/5 1/5 3/5 -3/5 -1/5 1/5 3/5 -3/5 -1/5  1/5  3/5 -3/5 -1/5  1/5  3/5 -3/5 -1/5 1/5 3/5 -3/5 -1/5 1/5 3/5;
                        -3/5 -3/5 -3/5 -3/5 -1/5 -1/5 -1/5 -1/5  1/5  1/5  1/5  1/5  3/5  3/5  3/5  3/5 -3/5 -3/5 -3/5 -3/5 -1/5 -1/5 -1/5 -1/5  1/5  1/5  1/5  1/5  3/5  3/5  3/5  3/5 -3/5 -3/5 -3/5 -3/5 -1/5 -1/5 -1/5 -1/5  1/5  1/5 1/5 1/5  3/5  3/5 3/5 3/5 -3/5 -3/5 -3/5 -3/5 -1/5 -1/5 -1/5 -1/5  1/5  1/5 1/5 1/5  3/5  3/5 3/5 3/5;
                        -3/5 -3/5 -3/5 -3/5 -3/5 -3/5 -3/5 -3/5 -3/5 -3/5 -3/5 -3/5 -3/5 -3/5 -3/5 -3/5 -1/5 -1/5 -1/5 -1/5 -1/5 -1/5 -1/5 -1/5 -1/5 -1/5 -1/5 -1/5 -1/5 -1/5 -1/5 -1/5  1/5  1/5  1/5  1/5  1/5  1/5  1/5  1/5  1/5  1/5 1/5 1/5  1/5  1/5 1/5 1/5  3/5  3/5  3/5  3/5  3/5  3/5  3/5  3/5  3/5  3/5 3/5 3/5  3/5  3/5 3/5 3/5];
               elseif( rem(n_ldof_ij^(1/3),1)==0 )
                 xi1 = linspace( -1, 1, n_ldof_ij^(1/3)+2 );
                 [xi2,xi1,xi3] = meshgrid( xi1(2:end-1), xi1(2:end-1), xi1(2:end-1) );
                 xii = [ xi1(:)'; xi2(:)'; xi3(:)' ];
               else
                 error(['lcoord_cell: hexahedral cells with ',num2str(n_ldof_ij),' internal dofs not supported.'])
               end

             end

         end

         xi = [ xi xii ];

     end

   end
 end