0001 function pe=PenaltyCreatePVector(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 pe=PenaltyCreatePVectorThermal(pengrid,grids,materials,inputs,kmax);
0014
0015 elseif strcmpi(analysis_type,'diagnostic_stokes')
0016
0017 pe=PenaltyCreatePVectorDiagnosticStokes(pengrid,grids,materials,inputs,kmax);
0018
0019 elseif strcmpi(analysis_type,'melting')
0020
0021 pe=PenaltyCreatePVectorMelting(pengrid,grids,materials,inputs,kmax);
0022
0023 end
0024 end
0025
0026 function pe=PenaltyCreatePVectorDiagnosticStokes(pengrid,grids,materials,inputs,kmax);
0027
0028 pe=elemvector(0);
0029
0030 end
0031
0032 function pe=PenaltyCreatePVectorThermal(pengrid,grids,materials,inputs,kmax);
0033
0034 if ~pengrid.active,
0035 pe=elemvector(0);
0036 else
0037
0038
0039 pe=elemvector(1);
0040
0041
0042 [pressure_param pressure_is_present]=recover_input(inputs,'pressure');
0043
0044 if (~pressure_is_present),
0045 error('CreatePVectorThermal error message: missing pressure as input');
0046 end
0047
0048
0049 doflist=grids(pengrid.grid).grid.doflist;
0050 dof=doflist(1);
0051 pe.row_indices=dof;
0052 pressure=pressure_param(dof);
0053
0054
0055 matpar=materials(end).constants;
0056 meltingpoint=matpar.meltingpoint;
0057 beta=matpar.beta;
0058 t_pmp=meltingpoint-beta*pressure;
0059
0060
0061 pe.terms(1)=kmax*10^(pengrid.penalty_offset)*t_pmp;
0062 end
0063 end
0064
0065
0066 function pe=PenaltyCreatePVectorMelting(pengrid,grids,materials,inputs,kmax);
0067
0068
0069 pe=elemvector(1);
0070
0071
0072 [pressure_param pressure_is_present]=recover_input(inputs,'pressure');
0073 [temperature_param temperature_is_present]=recover_input(inputs,'temperature');
0074 [melting_offset melting_offset_is_present]=recover_input(inputs,'melting_offset');
0075 if ~pengrid.thermal_steadystate,
0076 [dt dt_is_present]=recover_input(inputs,'dt');
0077 if ~dt_is_present
0078 error('CreatePVectorMelting error message: missing dt');
0079 end
0080 end
0081
0082 if (~temperature_is_present | ~pressure_is_present | ~melting_offset_is_present),
0083 error('CreatePVectorMelting error message: missing input parameters');
0084 end
0085
0086
0087 doflist=grids(pengrid.grid).grid.doflist;
0088 dof=doflist(1);
0089 pe.row_indices=dof;
0090 pressure=pressure_param(dof);
0091 temperature=temperature_param(dof);
0092
0093
0094 matpar=materials(end).constants;
0095 meltingpoint=matpar.meltingpoint;
0096 beta=matpar.beta;
0097 heatcapacity=matpar.heatcapacity;
0098 latentheat=matpar.latentheat;
0099
0100
0101 t_pmp=meltingpoint-beta*pressure;
0102
0103
0104
0105
0106
0107 if temperature<t_pmp
0108 pe.terms(1)=0;
0109 else
0110 if pengrid.thermal_steadystate,
0111 pe.terms(1)=melting_offset*10^(pengrid.penalty_offset)*(temperature-t_pmp);
0112 else
0113 pe.terms(1)=melting_offset*10^(pengrid.penalty_offset)*(temperature-t_pmp)/dt;
0114 end
0115 end
0116 end