0001 function Ke=PenaltyCreateKMatrix(pengrid,grids,materials,inputs,analysis_type,kmax)
0002
0003
0004
0005
0006
0007
0008
0009
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
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
0034 Ke=elemmatrix(1);
0035
0036
0037 doflist=grids(pengrid.grid).grid.doflist;
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
0045
0046 function Ke=CreateKMatrixMelting(pengrid,grids,materials,inputs,kmax)
0047
0048
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
0057 doflist=grids(pengrid.grid).grid.doflist;
0058 dof=doflist(1);
0059 pressure=pressure_param(dof);
0060 temperature=temperature_param(dof);
0061
0062
0063 matpar=materials(end).constants;
0064 meltingpoint=matpar.meltingpoint;
0065 beta=matpar.beta;
0066 t_pmp=meltingpoint-beta*pressure;
0067
0068
0069 if temperature>t_pmp
0070 Ke=elemmatrix(0);
0071
0072 else
0073
0074
0075 Ke=elemmatrix(1);
0076 Ke.row_indices=dof;
0077 Ke.terms(1,1)=kmax*10^pengrid.penalty_offset;
0078 end
0079
0080 end
0081
0082 function Ke=CreateKMatrixStokes(pengrid,grids,materials,inputs,kmax)
0083
0084
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
0092 slope=zeros(1,2);
0093 doflist=grids(pengrid.grid).grid.doflist;
0094 slope(1,1)=slope_param(doflist(1));
0095 slope(1,2)=slope_param(doflist(2));
0096
0097
0098 Ke=elemmatrix(3);
0099 for i=1:3,
0100 Ke.row_indices(i)=doflist(i);
0101 end
0102
0103
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