ModelProcessorMelting

PURPOSE ^

MODELPROCESSORMELTING - process model for a melting solution

SYNOPSIS ^

function [elements,grids,loads,constraints,materials,part,tpart]=ModelProcessorMelting(md,solutiontype);

DESCRIPTION ^

MODELPROCESSORMELTING - process model for a melting solution

   This routine uses all the informations in the model md and put them
   into different structures (grids, elements, loads, constrained,materials)
   that will be used to create the stiffness matrix and load vector.
   After this routine, the model md should not be called until the end of the
   solution sequence.

   Usage:
      [elements,grids,loads,constraints,materials,part,tpart]=ModelProcessorMelting(md,solutiontype)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [elements,grids,loads,constraints,materials,part,tpart]=ModelProcessorMelting(md,solutiontype);
0002 %MODELPROCESSORMELTING - process model for a melting solution
0003 %
0004 %   This routine uses all the informations in the model md and put them
0005 %   into different structures (grids, elements, loads, constrained,materials)
0006 %   that will be used to create the stiffness matrix and load vector.
0007 %   After this routine, the model md should not be called until the end of the
0008 %   solution sequence.
0009 %
0010 %   Usage:
0011 %      [elements,grids,loads,constraints,materials,part,tpart]=ModelProcessorMelting(md,solutiontype)
0012 
0013 global cluster
0014 
0015 if cluster,
0016     %We are running in parallel, we need to partition the elements
0017     element_partitioning=MeshPartition(md,numlabs);
0018 else
0019     %We are running in serial, all elements belong to the same partition.
0020     element_partitioning=ones(md.numberofelements,1);
0021     labindex=1; %older versions of matlab do not include the parallel toolbox labindex variable.
0022 end
0023 
0024 %Allocate grids and elements
0025 elements=struct('element',cell(md.numberofelements,1));
0026 materials=struct('material',cell(md.numberofelements+1,1));
0027 mygrids=zeros(md.numberofgrids,1); %this array determines grid partitioning.
0028 
0029 %Build elements
0030 if strcmpi(md.type,'2d'),
0031     error('ModelProcessorDiagnosticBaseVert error message: 2d mesh not supported yet!');
0032 end
0033 
0034 pos=find(element_partitioning==labindex);
0035 [elements(pos).element]=deal(pentaelem);
0036 
0037 elements(pos)=SetStructureField(elements(pos),'element','type','pentaelem');
0038 elements(pos)=SetStructureField(elements(pos),'element','id',pos);
0039 elements(pos)=SetStructureField(elements(pos),'element','g',md.elements(pos,1:6));
0040 elements(pos)=SetStructureField(elements(pos),'element','h',md.thickness(md.elements(pos,1:6)));
0041 elements(pos)=SetStructureField(elements(pos),'element','s',md.surface(md.elements(pos,1:6)));
0042 elements(pos)=SetStructureField(elements(pos),'element','b',md.bed(md.elements(pos,1:6)));
0043 elements(pos)=SetStructureField(elements(pos),'element','shelf',md.elementoniceshelf(pos));
0044 elements(pos)=SetStructureField(elements(pos),'element','onbed',md.elementonbed(pos));
0045 elements(pos)=SetStructureField(elements(pos),'element','onsurface',md.elementonsurface(pos));
0046 elements(pos)=SetStructureField(elements(pos),'element','acceleration',0);%acceleration 0 since no acceleration is posssible
0047 elements(pos)=SetStructureField(elements(pos),'element','collapse',1);
0048 elements(pos)=SetStructureField(elements(pos),'element','matid',pos);
0049 elements(pos)=SetStructureField(elements(pos),'element','melting',md.melting(md.elements(pos,1:6))/md.yts);
0050 elements(pos)=SetStructureField(elements(pos),'element','accumulation',md.accumulation(md.elements(pos,1:6))/md.yts);
0051 
0052 [materials(pos).material]=deal(matice);
0053 materials(pos)=SetStructureField(materials(pos),'material','id',pos);
0054 
0055 %Add physical constants in materials
0056 [materials(end).constants]=deal(matpar);
0057 materials(end)=SetStructureField(materials(end),'constants','g',md.g);
0058 materials(end)=SetStructureField(materials(end),'constants','rho_ice',md.rho_ice);
0059 materials(end)=SetStructureField(materials(end),'constants','rho_water',md.rho_water);
0060 materials(end)=SetStructureField(materials(end),'constants','thermalconductivity',md.thermalconductivity);
0061 materials(end)=SetStructureField(materials(end),'constants','heatcapacity',md.heatcapacity);
0062 materials(end)=SetStructureField(materials(end),'constants','latentheat',md.latentheat);
0063 materials(end)=SetStructureField(materials(end),'constants','beta',md.beta);
0064 materials(end)=SetStructureField(materials(end),'constants','meltingpoint',md.meltingpoint);
0065 materials(end)=SetStructureField(materials(end),'constants','mixed_layer_capacity',md.mixed_layer_capacity);
0066 materials(end)=SetStructureField(materials(end),'constants','thermal_exchange_velocity',md.thermal_exchange_velocity);
0067 
0068 if cluster, 
0069     %For elements, the corresponding grids belong to this cpu. Keep track of it.
0070     mygrids(md.elements(el3pos,:))=1;
0071     mygrids(md.elements(el6pos,:))=1;
0072     
0073     %Figure out which grids from the partitioning belong to different element partitions. We'll
0074     %call them 'border' grids.
0075     bordergrids=double(gplus(mygrids)>1);
0076 else
0077     bordergrids=zeros(md.numberofgrids,1); %no partitioning serially.
0078 end
0079 
0080 %Get the grids set up:
0081 grids=struct('grid',cell(md.numberofgrids,1));
0082 
0083 pos=[1:md.numberofgrids]';
0084 [grids(pos).grid]=deal(node);
0085 grids(pos)=SetStructureField(grids(pos),'grid','id',pos);
0086 grids(pos)=SetStructureField(grids(pos),'grid','x',md.x(pos));
0087 grids(pos)=SetStructureField(grids(pos),'grid','y',md.y(pos));
0088 grids(pos)=SetStructureField(grids(pos),'grid','z',md.z(pos));
0089 grids(pos)=SetStructureField(grids(pos),'grid','onbed',md.gridonbed(pos));
0090 grids(pos)=SetStructureField(grids(pos),'grid','border',bordergrids(pos));
0091 
0092 %spc degrees of freedom:
0093 %     for each grid, 6 degrees of freedom are allowed. These dofs are numbered from 1 to 6. The first 3
0094 %    deal with the (x,y,z) velocities, or deformations. The last 3 deal with the (x,y,z) rotations.
0095 %    If a certain degree of freedom i (1<=i<=6) is constrained to the value 0, the number i should be added to the
0096 %    gridset field of a grid.
0097 %    The gridset field holds all the numbers corresponding to the dofs that have been constrained to 0 value. Because
0098 %    we do not know firshand how many dofs have been constrained for a certain grid, we need a flexible way
0099 %    of keeping track of these constraints. Hence gridset is a string array, of no given size, with no given
0100 %    numerical order.
0101 %    Ex: if a grid has 0 values for the x and z deformations, and 0 values for the y rotation, we could add any of the
0102 %    following strings to the gridset: '135', '153', '315', etc ...
0103 grids(pos)=SetStructureField(grids(pos),'grid','gridset','123456');
0104 
0105 %for bed elements, free up the 3'rd dof for the first 3 grids
0106 for n=1:length(elements),
0107     %all dofs associated to this element are frozen, except if the element is on the bed, and acts as a 'fake' pentaelem,
0108     %and a true 'macayeal' tria element.
0109     if elements(n).element.onbed,
0110         for j=1:3,
0111             grids(elements(n).element.g(j)).grid.gridset='23456';
0112         end
0113     end
0114 end
0115 
0116 loads=struct('load',cell([0,1]));
0117 
0118 %create melting  penalties for every grid. The right hand side is computed using penalties
0119 count=1;
0120 for i=1:md.numberofgrids,
0121     if grids(i).grid.onbed,
0122         pengridobject=pengrid;
0123         pengridobject.id=count;
0124         pengridobject.dof=1;
0125         pengridobject.grid=i;
0126         if strcmpi(solutiontype,'meltingsteady'),
0127             pengridobject.thermal_steadystate=1;
0128         else
0129             pengridobject.thermal_steadystate=0;
0130         end
0131         pengridobject.penalty_offset=md.penalty_offset;
0132 
0133         loads(count).load=pengridobject;
0134         count=count+1;
0135     end
0136 end
0137 
0138 
0139 %Single point constraints: none;
0140 constraints=struct('constraint',cell(0,0));
0141 
0142 %Last thing, return a partitioning vector, and its transpose.
0143 [part,tpart]=PartitioningVector(md,grids);
0144 
0145 end %end function

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