PenaltyCreatePVector

PURPOSE ^

PENALTYCREATEPVECTOR - creates a penalty vector for a particular grid

SYNOPSIS ^

function pe=PenaltyCreatePVector(pengrid,grids,materials,inputs,analysis_type,kmax);

DESCRIPTION ^

PENALTYCREATEPVECTOR - creates a penalty vector for a particular grid

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

   Usage:
      pe=PenaltyCreatePVector(pengrid,grids,materials,inputs,analysis_type,kmax)
 
   See also CREATEPVECTOR, PENALTYCREATEKMATRIX, PENALTYCONTRAIN

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function pe=PenaltyCreatePVector(pengrid,grids,materials,inputs,analysis_type,kmax);
0002 %PENALTYCREATEPVECTOR - creates a penalty vector for a particular grid
0003 %
0004 %   This routine applies for thermal, melting, and Stokes solution.
0005 %
0006 %   Usage:
0007 %      pe=PenaltyCreatePVector(pengrid,grids,materials,inputs,analysis_type,kmax)
0008 %
0009 %   See also CREATEPVECTOR, PENALTYCREATEKMATRIX, PENALTYCONTRAIN
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 %end function
0025 
0026 function pe=PenaltyCreatePVectorDiagnosticStokes(pengrid,grids,materials,inputs,kmax);
0027     
0028     pe=elemvector(0);
0029 
0030 end %end function
0031 
0032 function pe=PenaltyCreatePVectorThermal(pengrid,grids,materials,inputs,kmax);
0033     
0034     if ~pengrid.active,
0035         pe=elemvector(0);
0036     else
0037         
0038         %initialize load vector
0039         pe=elemvector(1);
0040 
0041         %First recover pressure
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         %Recover pressure  at the grid and initialize elementary vector dof.
0049         doflist=grids(pengrid.grid).grid.doflist; %list of dofs in the g-set
0050         dof=doflist(1);
0051         pe.row_indices=dof;
0052         pressure=pressure_param(dof);
0053 
0054         %Compute pressure melting point
0055         matpar=materials(end).constants;
0056         meltingpoint=matpar.meltingpoint;
0057         beta=matpar.beta;
0058         t_pmp=meltingpoint-beta*pressure;
0059 
0060         %Add penalty load
0061         pe.terms(1)=kmax*10^(pengrid.penalty_offset)*t_pmp;
0062     end
0063 end %end function
0064 
0065 
0066 function pe=PenaltyCreatePVectorMelting(pengrid,grids,materials,inputs,kmax);
0067             
0068     %initialize load vector
0069     pe=elemvector(1);
0070 
0071     %First recover beta, pressure and temperature vectors;
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     %Recover pressure  at the grid and initialize elementary vector dof.
0087     doflist=grids(pengrid.grid).grid.doflist; %list of dofs in the g-set
0088     dof=doflist(1);
0089     pe.row_indices=dof;
0090     pressure=pressure_param(dof);
0091     temperature=temperature_param(dof);
0092 
0093     %Recover parameters
0094     matpar=materials(end).constants;
0095     meltingpoint=matpar.meltingpoint;
0096     beta=matpar.beta;
0097     heatcapacity=matpar.heatcapacity;
0098     latentheat=matpar.latentheat;
0099 
0100     %Compute pressure melting point
0101     t_pmp=meltingpoint-beta*pressure;
0102 
0103     %Add penalty load
0104     %This time, the penalty must have the same value as the one used for the thermal computation
0105     %so that the corresponding melting can be computed correctly
0106     %In the thermal computation, we used kmax=melting_offset, and the same penalty_offset
0107     if temperature<t_pmp %no melting
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 %end function

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