ModelProcessorDiagnosticStokes

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 numberofstokeselements=size(find(md.elements_type(:,2)==stokesenum()),1);
0026 elements=struct('element',cell(numberofstokeselements,1));
0027 materials=struct('material',cell(numberofstokeselements+1,1));
0028 mygrids=zeros(md.numberofgrids,1); %this array determines grid partitioning.
0029 
0030 %3d elements
0031 %pos=find(element_partitioning==labindex);
0032 pos=find(md.elements_type(:,2)==stokesenum());
0033 stokesnewnumber=[1:numberofstokeselements]';
0034 [elements(stokesnewnumber).element]=deal(pentaelem);
0035 
0036 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','type','pentaelem');
0037 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','id',stokesnewnumber);
0038 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','g',md.elements(pos,1:6));
0039 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','h',md.thickness(md.elements(pos,1:6)));
0040 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','s',md.surface(md.elements(pos,1:6)));
0041 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','b',md.bed(md.elements(pos,1:6)));
0042 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','friction_type',md.drag_type);
0043 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','k',md.drag(md.elements(pos,1:6)));
0044 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','p',md.p(pos));
0045 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','q',md.q(pos));
0046 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','shelf',md.elementoniceshelf(pos));
0047 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','onbed',md.elementonbed(pos));
0048 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','onsurface',md.elementonsurface(pos));
0049 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','meanvel',md.meanvel);
0050 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','epsvel',md.epsvel);
0051 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','acceleration',0);%acceleration 0 since no acceleration is posssible
0052 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','collapse',0);
0053 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','matid',stokesnewnumber);
0054 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','reconditioning_number',md.stokesreconditioning);
0055 
0056 [materials(stokesnewnumber).material]=deal(matice);
0057 
0058 materials(stokesnewnumber)=SetStructureField(materials(stokesnewnumber),'material','id',stokesnewnumber);
0059 materials(stokesnewnumber)=SetStructureField(materials(stokesnewnumber),'material','B',md.B(md.elements(pos,1:6))*[1;1;1;1;1;1]/6);
0060 materials(stokesnewnumber)=SetStructureField(materials(stokesnewnumber),'material','n',md.n(pos));
0061 
0062 
0063 %Add physical constants in materials
0064 [materials(end).constants]=deal(matpar);
0065 materials(end)=SetStructureField(materials(end),'constants','g',md.g);
0066 materials(end)=SetStructureField(materials(end),'constants','viscosity_overshoot',md.viscosity_overshoot);
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 
0078 if cluster, 
0079     %For elements, the corresponding grids belong to this cpu. Keep track of it.
0080     mygrids(md.elements(pos,:))=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','56');
0115 
0116 %Build grids on the border of stokes
0117 %Get border grids
0118 stokes_elements=find(md.elements_type(:,2)==stokesenum()); %find the elements on the stokes domain
0119 borderflags=zeros(md.numberofgrids,1); 
0120 borderflags(md.elements(stokes_elements,:))=1; %find all the grids of the elements on stokes domain, ie stokes grids and borderstokes
0121 borderstokes=borderflags-md.gridonstokes; %remove stokes grids from this list
0122 
0123 %spc the non stokes grids to 0 for all the dof
0124 for n=1:md.numberofgrids,
0125     if md.gridonstokes(n)==0, %the grids does not belong to the stokes part of the mesh
0126         grids(n).grid.gridset='123456'; %the six dofs are frozen
0127     end
0128     if borderstokes(n),
0129         grids(n).grid.gridset='12356';
0130     end
0131 end
0132 
0133 %Boundary conditions:
0134 
0135 %icefront
0136 segmentonneumann_diag_stokes=md.segmentonneumann_diag(find(md.elements_type(md.segmentonneumann_diag(:,end),2)),:);
0137 length_segmentonneumann_diag_stokes= size(segmentonneumann_diag_stokes,1);
0138 loads=struct('load',cell(length_segmentonneumann_diag_stokes,1));
0139 
0140 for i=1:length_segmentonneumann_diag_stokes,
0141 
0142     if (element_partitioning(segmentonneumann_diag_stokes(i,end))~=labindex), %this load does not belong to this cpu element partition.
0143         continue;
0144     end
0145     
0146     %build a quad ice front for the penta element
0147     loads(i).load=icefront;
0148     loads(i).load.eid=find(segmentonneumann_diag_stokes(i,end)==find(md.elements_type(:,2)==stokesenum())); %elements contain only stokes elements, so we have to renumbered
0149     loads(i).load.g(1)=segmentonneumann_diag_stokes(i,1);
0150     loads(i).load.g(2)=segmentonneumann_diag_stokes(i,2);
0151     loads(i).load.g(3)=segmentonneumann_diag_stokes(i,3);
0152     loads(i).load.g(4)=segmentonneumann_diag_stokes(i,4);
0153     loads(i).load.rho_water=md.rho_water;
0154     loads(i).load.type='quad';
0155 end
0156 
0157 %create penalties for grids on the base of icesheet. We must have wb=ub*db/dx+vb*db/dy
0158 count=length(loads)+1;
0159 for i=1:md.numberofgrids,
0160     if grids(i).grid.onbed & md.gridonicesheet(i)  & md.gridonstokes(i),
0161         pengridobject=pengrid;
0162         pengridobject.id=count;
0163         pengridobject.dof=1;
0164         pengridobject.grid=i;
0165         pengridobject.penalty_offset=md.penalty_offset;
0166 
0167         loads(count).load=pengridobject;
0168         count=count+1;
0169     end
0170 end
0171 
0172 %Single point constraints:
0173 spcs=find(md.gridondirichlet_diag);
0174 constraints=struct('constraint',cell(3*length(spcs),1));
0175 
0176 count=1;
0177 for i=1:md.numberofgrids,
0178     if ~md.gridonstokes(i),
0179 
0180         %constrain first dof
0181         constraints(count).constraint=spc;
0182         constraints(count).constraint.grid=i;
0183         constraints(count).constraint.dof=1;
0184         constraints(count).constraint.value=0; %this will be change to vx in the solution sequences
0185         count=count+1;
0186 
0187         %constrain second dof
0188         constraints(count).constraint=spc;
0189         constraints(count).constraint.grid=i;
0190         constraints(count).constraint.dof=2;
0191         constraints(count).constraint.value=0; %this will be change to vy in the solution sequences
0192         count=count+1;
0193 
0194         %constrain third dof
0195         constraints(count).constraint=spc;
0196         constraints(count).constraint.grid=i;
0197         constraints(count).constraint.dof=3;
0198         constraints(count).constraint.value=0; %this will be change to vz in the solution sequences
0199         count=count+1;
0200     end
0201 end
0202 
0203 
0204 %Last thing, return a partitioning vector, and its transpose.
0205 [part,tpart]=PartitioningVector(md,grids);
0206 
0207 
0208 end %end function

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