PenaltyCreateKMatrix

PURPOSE ^

PENALTYCREATEKMATRIX - creates a penalty matrix for a particular grid

SYNOPSIS ^

function Ke=PenaltyCreateKMatrix(pengrid,grids,materials,inputs,analysis_type,kmax)

DESCRIPTION ^

PENALTYCREATEKMATRIX - creates a penalty matrix for a particular grid

   This routine applies for thermal, melting, and Stokes solution.

   Usage:
      Ke=PenaltyCreateKMatrix(pengrid,grids,materials,inputs,analysis_type,kmax)
 
   See also CREATEPVECTOR, PENALTYCREATEPVECTOR, PENALTYCONTRAIN

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function Ke=PenaltyCreateKMatrix(pengrid,grids,materials,inputs,analysis_type,kmax)
0002 %PENALTYCREATEKMATRIX - creates a penalty matrix for a particular grid
0003 %
0004 %   This routine applies for thermal, melting, and Stokes solution.
0005 %
0006 %   Usage:
0007 %      Ke=PenaltyCreateKMatrix(pengrid,grids,materials,inputs,analysis_type,kmax)
0008 %
0009 %   See also CREATEPVECTOR, PENALTYCREATEPVECTOR, PENALTYCONTRAIN
0010 
0011     if (strcmpi(analysis_type,'thermaltransient') | strcmpi(analysis_type,'thermalsteady')),
0012 
0013         Ke=CreateKMatrixThermal(pengrid,grids,materials,inputs,kmax);
0014 
0015     elseif strcmpi(analysis_type,'diagnostic_stokes') 
0016 
0017         Ke=CreateKMatrixStokes(pengrid,grids,materials,inputs,kmax);
0018 
0019     elseif strcmpi(analysis_type,'melting')
0020 
0021         Ke=CreateKMatrixMelting(pengrid,grids,materials,inputs,kmax);
0022 
0023     end
0024 end %end function
0025 
0026 function Ke=CreateKMatrixThermal(pengrid,grids,materials,inputs,kmax)
0027 
0028     if ~pengrid.active,
0029         Ke=elemmatrix(0);
0030         return;
0031     else
0032         
0033         %initialize stiffness matrix
0034         Ke=elemmatrix(1);
0035 
0036         %Build linear indices for elementary stiffness matrix.
0037         doflist=grids(pengrid.grid).grid.doflist; %list of dofs in the g-set
0038         dof=doflist(1);
0039         Ke.row_indices=dof;
0040             
0041         Ke.terms(1,1)=kmax*10^(pengrid.penalty_offset);
0042     end
0043 
0044 end %function
0045 
0046 function Ke=CreateKMatrixMelting(pengrid,grids,materials,inputs,kmax)
0047     
0048     %First recover beta, pressure and temperature vectors;
0049     [pressure_param pressure_is_present]=recover_input(inputs,'pressure');
0050     [temperature_param temperature_is_present]=recover_input(inputs,'temperature');
0051 
0052     if (~temperature_is_present | ~pressure_is_present),
0053         error('CreatePVectorMelting error message: missing input parameters');
0054     end
0055 
0056     %Recover pressure  at the grid and initialize elementary vector dof.
0057     doflist=grids(pengrid.grid).grid.doflist; %list of dofs in the g-set
0058     dof=doflist(1);
0059     pressure=pressure_param(dof);
0060     temperature=temperature_param(dof);
0061 
0062     %Compute pressure melting point
0063     matpar=materials(end).constants;
0064     meltingpoint=matpar.meltingpoint;
0065     beta=matpar.beta;
0066     t_pmp=meltingpoint-beta*pressure;
0067 
0068     %Add penalty load
0069     if temperature>t_pmp %melting
0070         Ke=elemmatrix(0);
0071 
0072     else                 %no melting
0073         %If T<Tpmp, there must be no melting. Therefore, melting should be  constrained to 0 when T<Tpmp
0074         %instead of using spcs, one uses penalties
0075         Ke=elemmatrix(1);
0076         Ke.row_indices=dof;
0077         Ke.terms(1,1)=kmax*10^pengrid.penalty_offset;
0078     end    
0079 
0080 end %function
0081 
0082 function Ke=CreateKMatrixStokes(pengrid,grids,materials,inputs,kmax)
0083 
0084     %First recover beta, pressure and temperature vectors;
0085     [slope_param slope_is_present]=recover_input(inputs,'slopebed');
0086     
0087     if ~slope_is_present,
0088         error('CreatePVectorStokes error message: missing input parameters');
0089     end
0090     
0091     %Recover slope at the grid and initialize elementary vector dof.
0092     slope=zeros(1,2);
0093     doflist=grids(pengrid.grid).grid.doflist; %list of dofs in the g-set
0094     slope(1,1)=slope_param(doflist(1));
0095     slope(1,2)=slope_param(doflist(2));
0096 
0097     %Create elementary matrix
0098     Ke=elemmatrix(3);
0099     for i=1:3,
0100         Ke.row_indices(i)=doflist(i);
0101     end
0102 
0103     %Add penalty to contrain wb (wb=ub*db/dx+vb*db/dy)
0104     Ke.terms(3,1)=-slope(1,1)*kmax*10^pengrid.penalty_offset;
0105     Ke.terms(3,2)=-slope(1,2)*kmax*10^pengrid.penalty_offset;
0106     Ke.terms(3,3)=1*kmax*10^pengrid.penalty_offset;
0107 
0108 end %function

Generated on Sun 29-Mar-2009 20:22:55 by m2html © 2003