ModelProcessorDiagnosticHoriz

PURPOSE ^

MODELPROCESSORDIAGNOSTICHORIZ - process model for an horizontal diagnostic solution

SYNOPSIS ^

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

DESCRIPTION ^

MODELPROCESSORDIAGNOSTICHORIZ - process model for an horizontal diagnostic 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]=ModelProcessorDiagnosticHoriz(md)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [elements,grids,loads,constraints,materials,part,tpart]=ModelProcessorDiagnosticHoriz(md);
0002 %MODELPROCESSORDIAGNOSTICHORIZ - process model for an horizontal diagnostic 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]=ModelProcessorDiagnosticHoriz(md)
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 if strcmpi(md.type,'2d'),
0026     numberofelements_horiz=length(find(md.elements_type(:,1)==macayealenum()));
0027 else
0028     numberofelements_horiz=length(find(md.elements_type(:,1)==macayealenum() | md.elements_type(:,1)==pattynenum()));
0029 end
0030 elements=struct('element',cell(numberofelements_horiz,1));
0031 materials=struct('material',cell(numberofelements_horiz+1,1));
0032 mygrids=zeros(md.numberofgrids,1); %this array determines grid partitioning.
0033 
0034 %Deal with 2d elements
0035 if strcmpi(md.type,'2d'),
0036     el3pos=find((md.elements_type(:,1)==macayealenum()) & (element_partitioning==labindex));
0037     if ~isempty(el3pos),
0038         pos=[1:length(el3pos)]';
0039         [elements(pos).element]=deal(triaelem);
0040 
0041         elements(pos)=SetStructureField(elements(pos),'element','type','triaelem');
0042         elements(pos)=SetStructureField(elements(pos),'element','id',pos);
0043         elements(pos)=SetStructureField(elements(pos),'element','g',md.elements(el3pos,:));
0044         elements(pos)=SetStructureField(elements(pos),'element','h',md.thickness(md.elements(el3pos,1:3)));
0045         elements(pos)=SetStructureField(elements(pos),'element','s',md.surface(md.elements(el3pos,1:3)));
0046         elements(pos)=SetStructureField(elements(pos),'element','b',md.bed(md.elements(el3pos,1:3)));
0047         elements(pos)=SetStructureField(elements(pos),'element','friction_type',md.drag_type);
0048         elements(pos)=SetStructureField(elements(pos),'element','k',md.drag(md.elements(el3pos,1:3)));
0049         elements(pos)=SetStructureField(elements(pos),'element','p',md.p(el3pos));
0050         elements(pos)=SetStructureField(elements(pos),'element','q',md.q(el3pos));
0051         elements(pos)=SetStructureField(elements(pos),'element','shelf',md.elementoniceshelf(el3pos));
0052         elements(pos)=SetStructureField(elements(pos),'element','meanvel',md.meanvel);
0053         elements(pos)=SetStructureField(elements(pos),'element','epsvel',md.epsvel);
0054         elements(pos)=SetStructureField(elements(pos),'element','acceleration',md.acceleration);
0055         elements(pos)=SetStructureField(elements(pos),'element','matid',pos);
0056 
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(el3pos,1:3))*[1;1;1]/3);
0062         materials(pos)=SetStructureField(materials(pos),'material','n',md.n(el3pos));
0063     end
0064 else
0065 
0066     %3d elements
0067     %First create penta elements (ie non hutter)
0068     el6pos=find((md.elements_type(:,1)==macayealenum() | md.elements_type(:,1)==pattynenum()) & (element_partitioning==labindex));
0069 
0070     if ~isempty(el6pos),
0071         pos=[1:length(el6pos)]'; %Renumber the elements
0072         [elements(pos).element]=deal(pentaelem);
0073 
0074         elements(pos)=SetStructureField(elements(pos),'element','type','pentaelem');
0075         elements(pos)=SetStructureField(elements(pos),'element','id',pos);
0076         elements(pos)=SetStructureField(elements(pos),'element','g',md.elements(el6pos,1:6));
0077         elements(pos)=SetStructureField(elements(pos),'element','h',md.thickness(md.elements(el6pos,1:6)));
0078         elements(pos)=SetStructureField(elements(pos),'element','s',md.surface(md.elements(el6pos,1:6)));
0079         elements(pos)=SetStructureField(elements(pos),'element','b',md.bed(md.elements(el6pos,1:6)));
0080         elements(pos)=SetStructureField(elements(pos),'element','friction_type',md.drag_type);
0081         elements(pos)=SetStructureField(elements(pos),'element','k',md.drag(md.elements(el6pos,1:6)));
0082         elements(pos)=SetStructureField(elements(pos),'element','p',md.p(el6pos));
0083         elements(pos)=SetStructureField(elements(pos),'element','q',md.q(el6pos));
0084         elements(pos)=SetStructureField(elements(pos),'element','shelf',md.elementoniceshelf(el6pos));
0085         elements(pos)=SetStructureField(elements(pos),'element','onbed',md.elementonbed(el6pos));
0086         elements(pos)=SetStructureField(elements(pos),'element','onsurface',md.elementonsurface(el6pos));
0087         elements(pos)=SetStructureField(elements(pos),'element','meanvel',md.meanvel);
0088         elements(pos)=SetStructureField(elements(pos),'element','epsvel',md.epsvel);
0089         elements(pos)=SetStructureField(elements(pos),'element','acceleration',0);%acceleration 0 since no acceleration is posssible
0090         elements(pos)=SetStructureField(elements(pos),'element','collapse',0);
0091         elements(pos)=SetStructureField(elements(pos),'element','matid',pos);
0092 
0093         [materials(pos).material]=deal(matice);
0094 
0095         materials(pos)=SetStructureField(materials(pos),'material','id',pos);
0096         materials(pos)=SetStructureField(materials(pos),'material','B',md.B(md.elements(el6pos,1:6))*[1;1;1;1;1;1]/6);
0097         materials(pos)=SetStructureField(materials(pos),'material','n',md.n(el6pos));
0098 
0099         %For penta elements where we want to implement MacAyeal's element, we need to collapse
0100         %the formulation into trias:
0101         pos=find(ismember(el6pos,find(md.elements_type(:,1)==macayealenum())));
0102         elements(pos)=SetStructureField(elements(pos),'element','collapse',ones(length(pos),1));
0103 
0104     end
0105 
0106 end
0107 
0108 
0109 %Add physical constants in materials
0110 [materials(end).constants]=deal(matpar);
0111 materials(end)=SetStructureField(materials(end),'constants','g',md.g);
0112 materials(end)=SetStructureField(materials(end),'constants','viscosity_overshoot',md.viscosity_overshoot);
0113 materials(end)=SetStructureField(materials(end),'constants','rho_ice',md.rho_ice);
0114 materials(end)=SetStructureField(materials(end),'constants','rho_water',md.rho_water);
0115 materials(end)=SetStructureField(materials(end),'constants','thermalconductivity',md.thermalconductivity);
0116 materials(end)=SetStructureField(materials(end),'constants','heatcapacity',md.heatcapacity);
0117 materials(end)=SetStructureField(materials(end),'constants','latentheat',md.latentheat);
0118 materials(end)=SetStructureField(materials(end),'constants','beta',md.beta);
0119 materials(end)=SetStructureField(materials(end),'constants','meltingpoint',md.meltingpoint);
0120 materials(end)=SetStructureField(materials(end),'constants','mixed_layer_capacity',md.mixed_layer_capacity);
0121 materials(end)=SetStructureField(materials(end),'constants','thermal_exchange_velocity',md.thermal_exchange_velocity);
0122 
0123 if cluster, 
0124     %For elements, the corresponding grids belong to this cpu. Keep track of it.
0125     mygrids(md.elements(el3pos,:))=1;
0126     mygrids(md.elements(el6pos,:))=1;
0127 end
0128 
0129 %Deal with acceleration MacAyeal's element.
0130 if md.acceleration & strcmpi(md.type,'2d'), 
0131     %build accelerated element.
0132     elements(end+1).element=acceleratedtriaelem;
0133     elements(end).element.type='acceleratedtriaelem';
0134     elements(end).element.x=md.x;
0135     elements(end).element.y=md.y;
0136     elements(end).element.index=md.elements;
0137     elements(end).element.nods=md.numberofgrids;
0138     elements(end).element.nel=md.numberofelements;
0139     elements(end).element.thickness=md.thickness;
0140     elements(end).element.surface=md.surface;
0141     elements(end).element.bed=md.bed;
0142     elements(end).element.thickness_el=(md.thickness(md.elements))*[1;1;1]/3;
0143     elements(end).element.B_bar=(md.B(md.elements))*[1;1;1]/3;
0144     elements(end).element.glen_coeff=md.n;
0145     elements(end).element.friction_type=md.drag_type;
0146     elements(end).element.drag=md.drag;
0147     elements(end).element.p=md.p;
0148     elements(end).element.q=md.q;
0149     elements(end).element.meanvel=md.meanvel;
0150     elements(end).element.epsvel=md.epsvel;
0151     if isempty(md.segmentonneumann_diag),
0152         elements(end).element.index_icefront=[];    
0153     else
0154         elements(end).element.index_icefront=md.segmentonneumann_diag(:,1:2);
0155     end
0156     elements(end).element.gridoniceshelf=md.gridoniceshelf;
0157     elements(end).element.elementonicesheet=md.elementonicesheet;
0158     elements(end).element.matid=1; %point to first matid element.
0159 end
0160 
0161 if cluster, 
0162     %Figure out which grids from the partitioning belong to different element partitions. We'll
0163     %call them 'border' grids.
0164     bordergrids=double(gplus(mygrids)>1);
0165 else
0166     bordergrids=zeros(md.numberofgrids,1); %no partitioning serially.
0167 end
0168 
0169 %Get the grids set up:
0170 grids=struct('grid',cell(md.numberofgrids,1));
0171 
0172 pos=[1:md.numberofgrids]';
0173 [grids(pos).grid]=deal(node);
0174 grids(pos)=SetStructureField(grids(pos),'grid','id',pos);
0175 grids(pos)=SetStructureField(grids(pos),'grid','x',md.x(pos));
0176 grids(pos)=SetStructureField(grids(pos),'grid','y',md.y(pos));
0177 grids(pos)=SetStructureField(grids(pos),'grid','z',md.z(pos));
0178 grids(pos)=SetStructureField(grids(pos),'grid','onbed',md.gridonbed(pos));
0179 grids(pos)=SetStructureField(grids(pos),'grid','border',bordergrids(pos));
0180 
0181 %spc degrees of freedom:
0182 %     for each grid, 6 degrees of freedom are allowed. These dofs are numbered from 1 to 6. The first 3
0183 %    deal with the (x,y,z) velocities, or deformations. The last 3 deal with the (x,y,z) rotations.
0184 %    If a certain degree of freedom i (1<=i<=6) is constrained to the value 0, the number i should be added to the
0185 %    gridset field of a grid.
0186 %    The gridset field holds all the numbers corresponding to the dofs that have been constrained to 0 value. Because
0187 %    we do not know firshand how many dofs have been constrained for a certain grid, we need a flexible way
0188 %    of keeping track of these constraints. Hence gridset is a string array, of no given size, with no given
0189 %    numerical order.
0190 %    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
0191 %    following strings to the gridset: '135', '153', '315', etc ...
0192 if strcmpi(md.type,'3d')
0193     grids(pos)=SetStructureField(grids(pos),'grid','gridset','123456');
0194 else
0195     grids(pos)=SetStructureField(grids(pos),'grid','gridset','3456');
0196 end
0197 
0198 %spc macayeal type grids  in 3d meshes
0199 if strcmpi(md.type,'3d'),
0200     for n=1:length(elements),
0201         if md.elements_type(el6pos(n),1)==macayealenum(),
0202             %all dofs associated to this element are frozen, except if the element is on the bed, and acts as a 'fake' pentaelem,
0203             %and a true 'macayeal' tria element.
0204             if elements(n).element.onbed,
0205                 for j=1:3,
0206                     grids(elements(n).element.g(j)).grid.gridset='3456';
0207                 end
0208             end
0209         else
0210             for j=1:length(elements(n).element.g),
0211                 grids(elements(n).element.g(j)).grid.gridset='3456';
0212             end
0213         end
0214     end
0215 end
0216 
0217 %Boundary conditions:
0218 
0219 %icefront
0220 loads=struct('load',cell(length(md.segmentonneumann_diag),1));
0221 
0222 for i=1:size(md.segmentonneumann_diag,1),
0223 
0224     if (element_partitioning(md.segmentonneumann_diag(i,end))~=labindex), %this load does not belong to this cpu element partition.
0225         continue;
0226     end
0227 
0228     if strcmpi(md.type,'3d'),
0229         if md.elements_type(md.segmentonneumann_diag(i,end))==macayealenum(),
0230             loads(i).load=icefront;
0231             loads(i).load.eid=find(el6pos==md.segmentonneumann_diag(i,end));
0232             loads(i).load.g(1)=md.segmentonneumann_diag(i,1);
0233             loads(i).load.g(2)=md.segmentonneumann_diag(i,2);
0234             loads(i).load.rho_water=md.rho_water;
0235             loads(i).load.type='segment';
0236 
0237         elseif md.elements_type(md.segmentonneumann_diag(i,end))==pattynenum(),
0238             %build a quad ice front for the penta element
0239             loads(i).load=icefront;
0240             loads(i).load.eid=find(el6pos==md.segmentonneumann_diag(i,end));
0241             loads(i).load.g(1)=md.segmentonneumann_diag(i,1);
0242             loads(i).load.g(2)=md.segmentonneumann_diag(i,2);
0243             loads(i).load.g(3)=md.segmentonneumann_diag(i,3);
0244             loads(i).load.g(4)=md.segmentonneumann_diag(i,4);
0245             loads(i).load.type='quad';
0246 
0247         else
0248             if ~(md.elements_type(md.segmentonneumann_diag(i,end))==hutterenum()),
0249                 error('diagnostic error message: unsupported  element type');
0250             end
0251         end
0252     else
0253         if md.elements_type(md.segmentonneumann_diag(i,end))==macayealenum(),
0254             loads(i).load=icefront;
0255             loads(i).load.eid=find(el3pos==md.segmentonneumann_diag(i,end));
0256             loads(i).load.g(1)=md.segmentonneumann_diag(i,1);
0257             loads(i).load.g(2)=md.segmentonneumann_diag(i,2);
0258             loads(i).load.rho_water=md.rho_water;
0259             loads(i).load.type='segment';
0260 
0261         else
0262             if ~(md.elements_type(md.segmentonneumann_diag(i,end))==hutterenum()),
0263                 error('diagnostic error message: unsupported  element type');
0264             end
0265         end
0266     
0267     end
0268 
0269 end
0270 
0271 %Initialize constraints structure
0272 constraints=struct('constraint',cell(2*length(find(md.gridondirichlet_diag)),1));
0273 
0274 count=1;
0275 for i=1:md.numberofgrids,
0276     %deal with spcs
0277     if md.gridonhutter(i),
0278 
0279         %constrain first dof
0280         constraints(count).constraint=spc;
0281         constraints(count).constraint.grid=i;
0282         constraints(count).constraint.dof=1;
0283         constraints(count).constraint.value=0; %this will be change to vx in the solution sequences
0284         count=count+1;
0285 
0286         %constrain second dof
0287         constraints(count).constraint=spc;
0288         constraints(count).constraint.grid=i;
0289         constraints(count).constraint.dof=2;
0290         constraints(count).constraint.value=0; %this will be change to vy in the solution sequences
0291         count=count+1;
0292 
0293     elseif (md.gridonmacayeal(i) | md.gridonpattyn(i)) &  md.gridondirichlet_diag(i),
0294 
0295         %constrain first dof
0296         constraints(count).constraint=spc;
0297         constraints(count).constraint.grid=i;
0298         constraints(count).constraint.dof=1;
0299         constraints(count).constraint.value=md.dirichletvalues_diag(i,1)/md.yts; %in m/s
0300         count=count+1;
0301 
0302         %constrain second dof
0303         constraints(count).constraint=spc;
0304         constraints(count).constraint.grid=i;
0305         constraints(count).constraint.dof=2;
0306         constraints(count).constraint.value=md.dirichletvalues_diag(i,2)/md.yts; %in m/s
0307         count=count+1;
0308     end
0309 end
0310 
0311 %deal with mpcs
0312 if ~isempty(md.penalties) & ~isnan(md.penalties),
0313     for i=1:size(md.penalties,1),
0314         for j=1:(md.numlayers-1),
0315 
0316             %constrain first dof
0317             constraints(count).constraint=rgb;
0318             constraints(count).constraint.grid1=md.penalties(i,1);
0319             constraints(count).constraint.grid2=md.penalties(i,j+1);
0320             constraints(count).constraint.dof=1;
0321             count=count+1;
0322             
0323             %constrain second dof
0324             constraints(count).constraint=rgb;
0325             constraints(count).constraint.grid1=md.penalties(i,1);
0326             constraints(count).constraint.grid2=md.penalties(i,j+1);
0327             constraints(count).constraint.dof=2;
0328             count=count+1;
0329         end
0330     end
0331 end
0332 
0333 %Last thing, return a partitioning vector, and its transpose.
0334 [part,tpart]=PartitioningVector(md,grids);
0335 
0336 end %end function

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