ModelProcessorThermal

PURPOSE ^

MODELPROCESSORTHERMAL - process model for a thermal solution

SYNOPSIS ^

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

DESCRIPTION ^

MODELPROCESSORTHERMAL - process model for a thermal 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]=ModelProcessorThermal(md,solutiontype)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [elements,grids,loads,constraints,materials,part,tpart]=ModelProcessorThermal(md,solutiontype);
0002 %MODELPROCESSORTHERMAL - process model for a thermal 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]=ModelProcessorThermal(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 3d elements
0030 pos=find(element_partitioning==labindex);
0031 [elements(pos).element]=deal(pentaelem);
0032 
0033 elements(pos)=SetStructureField(elements(pos),'element','type','pentaelem');
0034 elements(pos)=SetStructureField(elements(pos),'element','id',pos);
0035 elements(pos)=SetStructureField(elements(pos),'element','g',md.elements(pos,1:6));
0036 elements(pos)=SetStructureField(elements(pos),'element','h',md.thickness(md.elements(pos,1:6)));
0037 elements(pos)=SetStructureField(elements(pos),'element','s',md.surface(md.elements(pos,1:6)));
0038 elements(pos)=SetStructureField(elements(pos),'element','b',md.bed(md.elements(pos,1:6)));
0039 elements(pos)=SetStructureField(elements(pos),'element','friction_type',md.drag_type);
0040 elements(pos)=SetStructureField(elements(pos),'element','k',md.drag(md.elements(pos,1:6)));
0041 elements(pos)=SetStructureField(elements(pos),'element','p',md.p(pos));
0042 elements(pos)=SetStructureField(elements(pos),'element','q',md.q(pos));
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','meanvel',md.meanvel);
0047 elements(pos)=SetStructureField(elements(pos),'element','epsvel',md.epsvel);
0048 elements(pos)=SetStructureField(elements(pos),'element','acceleration',0);%acceleration 0 since no acceleration is posssible
0049 elements(pos)=SetStructureField(elements(pos),'element','collapse',0);
0050 elements(pos)=SetStructureField(elements(pos),'element','matid',pos);
0051 elements(pos)=SetStructureField(elements(pos),'element','geothermalflux',md.geothermalflux(md.elements(pos,1:6)));
0052 if strcmpi(solutiontype,'thermalsteady'),
0053     elements(pos)=SetStructureField(elements(pos),'element','thermal_steadystate',1);
0054 else
0055     elements(pos)=SetStructureField(elements(pos),'element','thermal_steadystate',0);
0056 end
0057 
0058 [materials(pos).material]=deal(matice);
0059 
0060 materials(pos)=SetStructureField(materials(pos),'material','id',pos);
0061 materials(pos)=SetStructureField(materials(pos),'material','B',md.B(md.elements(pos,1:6))*[1;1;1;1;1;1]/6);
0062 materials(pos)=SetStructureField(materials(pos),'material','n',md.n(pos));
0063 
0064 %Add physical constants in materials
0065 [materials(end).constants]=deal(matpar);
0066 materials(end)=SetStructureField(materials(end),'constants','g',md.g);
0067 materials(end)=SetStructureField(materials(end),'constants','rho_ice',md.rho_ice);
0068 materials(end)=SetStructureField(materials(end),'constants','rho_water',md.rho_water);
0069 materials(end)=SetStructureField(materials(end),'constants','thermalconductivity',md.thermalconductivity);
0070 materials(end)=SetStructureField(materials(end),'constants','heatcapacity',md.heatcapacity);
0071 materials(end)=SetStructureField(materials(end),'constants','latentheat',md.latentheat);
0072 materials(end)=SetStructureField(materials(end),'constants','beta',md.beta);
0073 materials(end)=SetStructureField(materials(end),'constants','meltingpoint',md.meltingpoint);
0074 materials(end)=SetStructureField(materials(end),'constants','mixed_layer_capacity',md.mixed_layer_capacity);
0075 materials(end)=SetStructureField(materials(end),'constants','thermal_exchange_velocity',md.thermal_exchange_velocity);
0076 
0077 if cluster, 
0078     %For elements, the corresponding grids belong to this cpu. Keep track of it.
0079     mygrids(md.elements(el3pos,:))=1;
0080     mygrids(md.elements(el6pos,:))=1;
0081 end
0082 
0083 if cluster, 
0084     %Figure out which grids from the partitioning belong to different element partitions. We'll
0085     %call them 'border' grids.
0086     bordergrids=double(gplus(mygrids)>1);
0087 else
0088     bordergrids=zeros(md.numberofgrids,1); %no partitioning serially.
0089 end
0090 
0091 %Get the grids set up:
0092 grids=struct('grid',cell(md.numberofgrids,1));
0093 
0094 pos=[1:md.numberofgrids]';
0095 [grids(pos).grid]=deal(node);
0096 grids(pos)=SetStructureField(grids(pos),'grid','id',pos);
0097 grids(pos)=SetStructureField(grids(pos),'grid','x',md.x(pos));
0098 grids(pos)=SetStructureField(grids(pos),'grid','y',md.y(pos));
0099 grids(pos)=SetStructureField(grids(pos),'grid','z',md.z(pos));
0100 grids(pos)=SetStructureField(grids(pos),'grid','onbed',md.gridonbed(pos));
0101 grids(pos)=SetStructureField(grids(pos),'grid','border',bordergrids(pos));
0102 
0103 %spc degrees of freedom:
0104 %     for each grid, 6 degrees of freedom are allowed. These dofs are numbered from 1 to 6. The first 3
0105 %    deal with the (x,y,z) velocities, or deformations. The last 3 deal with the (x,y,z) rotations.
0106 %    If a certain degree of freedom i (1<=i<=6) is constrained to the value 0, the number i should be added to the
0107 %    gridset field of a grid.
0108 %    The gridset field holds all the numbers corresponding to the dofs that have been constrained to 0 value. Because
0109 %    we do not know firshand how many dofs have been constrained for a certain grid, we need a flexible way
0110 %    of keeping track of these constraints. Hence gridset is a string array, of no given size, with no given
0111 %    numerical order.
0112 %    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
0113 %    following strings to the gridset: '135', '153', '315', etc ...
0114 grids(pos)=SetStructureField(grids(pos),'grid','gridset','23456');
0115 
0116 %Boundary conditions:
0117 
0118 %penalties
0119 loads=struct('load',cell(0,1));
0120 
0121 %create thermal penalties for every grid. Temperatures can never go over the melting point.
0122 count=1;
0123 for i=1:md.numberofgrids,
0124     if ~md.gridondirichlet_thermal(i),  %No penalty applied on spc grids!
0125         pengridobject=pengrid;
0126         pengridobject.id=count;
0127         pengridobject.dof=1;
0128         pengridobject.grid=i;
0129         pengridobject.active=0;
0130         if strcmpi(solutiontype,'thermalsteady'),
0131             pengridobject.thermal_steadystate=1;
0132         else
0133             pengridobject.thermal_steadystate=0;
0134         end
0135         pengridobject.penalty_offset=md.penalty_offset;
0136 
0137         loads(count).load=pengridobject;
0138         count=count+1;
0139     end
0140 end
0141 
0142 %Single point constraints:
0143 spcs=find(md.gridondirichlet_thermal);
0144 constraints=struct('constraint',cell(length(spcs),1));
0145 
0146 count=1;
0147 for i=1:md.numberofgrids,
0148     if md.gridondirichlet_thermal(i),
0149         constraints(count).constraint=spc;
0150         constraints(count).constraint.grid=i;
0151         constraints(count).constraint.dof=1;
0152         constraints(count).constraint.value=md.dirichletvalues_thermal(i,1);
0153         count=count+1;
0154     end
0155 end
0156 
0157 
0158 %Last thing, return a partitioning vector, and its transpose.
0159 [part,tpart]=PartitioningVector(md,grids);
0160 
0161 
0162 end %end function

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