Index: /issm/trunk/src/m/solutions/ice/ModelProcessor.m =================================================================== --- /issm/trunk/src/m/solutions/ice/ModelProcessor.m (revision 326) +++ /issm/trunk/src/m/solutions/ice/ModelProcessor.m (revision 327) @@ -18,6 +18,4 @@ elseif strcmpi(solutiontype,'diagnostic_hutter'), [elements,grids,loads,constraints,materials,part,tpart]=ModelProcessorDiagnosticHutter(md); -elseif strcmpi(solutiontype,'diagnostic_basevert'), - [elements,grids,loads,constraints,materials,part,tpart]=ModelProcessorDiagnosticBaseVert(md); elseif strcmpi(solutiontype,'thermalsteady')| strcmpi(solutiontype,'thermaltransient'), [elements,grids,loads,constraints,materials,part,tpart]=ModelProcessorThermal(md,solutiontype); Index: sm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticBaseVert.m =================================================================== --- /issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticBaseVert.m (revision 326) +++ (revision ) @@ -1,128 +1,0 @@ -function [elements,grids,loads,constraints,materials,part,tpart]=ModelProcessorDiagnosticBaseVert(md); -%MODELPROCESSORDIAGNOSTICBASEVERT - process the model for a DiagmosticBaseVert 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]=ModelProcessorDiagnosticBaseVert(md); - -global cluster - -if cluster, - %We are running in parallel, we need to partition the elements - element_partitioning=MeshPartition(md,numlabs); -else - %We are running in serial, all elements belong to the same partition. - element_partitioning=ones(md.numberofelements,1); - labindex=1; %older versions of matlab do not include the parallel toolbox labindex variable. -end - -%Allocate grids and elements -elements=struct('element',cell(md.numberofelements,1)); -materials=struct('material',cell(md.numberofelements+1,1)); -mygrids=zeros(md.numberofgrids,1); %this array determines grid partitioning. - -%Build elements -if strcmpi(md.type,'2d'), - error('ModelProcessorDiagnosticBaseVert error message: 2d mesh not supported yet!'); -end - -pos=find(element_partitioning==labindex); -[elements(pos).element]=deal(pentaelem); - -elements(pos)=SetStructureField(elements(pos),'element','type','pentaelem'); -elements(pos)=SetStructureField(elements(pos),'element','id',pos); -elements(pos)=SetStructureField(elements(pos),'element','g',md.elements(pos,1:6)); -elements(pos)=SetStructureField(elements(pos),'element','h',md.thickness(md.elements(pos,1:6))); -elements(pos)=SetStructureField(elements(pos),'element','s',md.surface(md.elements(pos,1:6))); -elements(pos)=SetStructureField(elements(pos),'element','b',md.bed(md.elements(pos,1:6))); -elements(pos)=SetStructureField(elements(pos),'element','shelf',md.elementoniceshelf(pos)); -elements(pos)=SetStructureField(elements(pos),'element','onbed',md.elementonbed(pos)); -elements(pos)=SetStructureField(elements(pos),'element','onsurface',md.elementonsurface(pos)); -elements(pos)=SetStructureField(elements(pos),'element','collapse',1); -elements(pos)=SetStructureField(elements(pos),'element','matid',pos); -elements(pos)=SetStructureField(elements(pos),'element','melting',full(md.melting(md.elements(pos,1:6))/md.yts)); -elements(pos)=SetStructureField(elements(pos),'element','accumulation',md.accumulation(md.elements(pos,1:6))/md.yts); - -[materials(pos).material]=deal(matice); - -materials(pos)=SetStructureField(materials(pos),'material','id',pos); - -%Add physical constants in materials -[materials(end).constants]=deal(matpar); -materials(end)=SetStructureField(materials(end),'constants','g',md.g); -materials(end)=SetStructureField(materials(end),'constants','rho_ice',md.rho_ice); -materials(end)=SetStructureField(materials(end),'constants','rho_water',md.rho_water); -materials(end)=SetStructureField(materials(end),'constants','thermalconductivity',md.thermalconductivity); -materials(end)=SetStructureField(materials(end),'constants','heatcapacity',md.heatcapacity); -materials(end)=SetStructureField(materials(end),'constants','latentheat',md.latentheat); -materials(end)=SetStructureField(materials(end),'constants','beta',md.beta); -materials(end)=SetStructureField(materials(end),'constants','meltingpoint',md.meltingpoint); -materials(end)=SetStructureField(materials(end),'constants','mixed_layer_capacity',md.mixed_layer_capacity); -materials(end)=SetStructureField(materials(end),'constants','thermal_exchange_velocity',md.thermal_exchange_velocity); - - -if cluster, - %For elements, the corresponding grids belong to this cpu. Keep track of it. - mygrids(md.elements(el3pos,:))=1; - mygrids(md.elements(el6pos,:))=1; - - %Figure out which grids from the partitioning belong to different element partitions. We'll - %call them 'border' grids. - bordergrids=double(gplus(mygrids)>1); -else - bordergrids=zeros(md.numberofgrids,1); %no partitioning serially. -end - -%Get the grids set up: -grids=struct('grid',cell(md.numberofgrids,1)); - -pos=[1:md.numberofgrids]'; -[grids(pos).grid]=deal(node); -grids(pos)=SetStructureField(grids(pos),'grid','id',pos); -grids(pos)=SetStructureField(grids(pos),'grid','x',md.x(pos)); -grids(pos)=SetStructureField(grids(pos),'grid','y',md.y(pos)); -grids(pos)=SetStructureField(grids(pos),'grid','z',md.z(pos)); -grids(pos)=SetStructureField(grids(pos),'grid','onbed',md.gridonbed(pos)); -grids(pos)=SetStructureField(grids(pos),'grid','border',bordergrids(pos)); - -%spc degrees of freedom: -% for each grid, 6 degrees of freedom are allowed. These dofs are numbered from 1 to 6. The first 3 -% deal with the (x,y,z) velocities, or deformations. The last 3 deal with the (x,y,z) rotations. -% If a certain degree of freedom i (1<=i<=6) is constrained to the value 0, the number i should be added to the -% gridset field of a grid. -% The gridset field holds all the numbers corresponding to the dofs that have been constrained to 0 value. Because -% we do not know firshand how many dofs have been constrained for a certain grid, we need a flexible way -% of keeping track of these constraints. Hence gridset is a string array, of no given size, with no given -% numerical order. -% 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 -% following strings to the gridset: '135', '153', '315', etc ... - -%Here, block all dof for every grid -grids(pos)=SetStructureField(grids(pos),'grid','gridset','123456'); - -%for bed elements, free up the 3'rd dof for the first 3 grids -for n=1:length(elements), - %all dofs associated to this element are frozen, except if the element is on the bed, and acts as a 'fake' pentaelem, - %and a true 'macayeal' tria element. - if elements(n).element.onbed, - for j=1:3, - grids(elements(n).element.g(j)).grid.gridset='12456'; - end - end -end - -%Boundary conditions: no penalties applied here -loads=struct('load',cell([0,1])); - -%Single point constraints: none; -constraints=struct('constraint',cell(0,0)); - -%Last thing, return a partitioning vector, and its transpose. -[part,tpart]=PartitioningVector(md,grids); - -end %end function Index: /issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticVert.m =================================================================== --- /issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticVert.m (revision 326) +++ /issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticVert.m (revision 327) @@ -108,20 +108,6 @@ loads=struct('load',cell([0,1])); -%Single point constraints: -spcs=find(md.gridonbed); -constraints=struct('constraint',cell(length(spcs),1)); - -count=1; -for i=1:md.numberofgrids, - if md.gridonbed(i), - - %constrain first dof - constraints(count).constraint=spc; - constraints(count).constraint.grid=i; - constraints(count).constraint.dof=3; - constraints(count).constraint.value=0; %this will be change to vz in the solution sequences. - count=count+1; - end -end +%No single constraints +constraints=struct('constraint',cell(0,0)); %Last thing, return a partitioning vector, and its transpose. Index: /issm/trunk/src/m/solutions/ice/diagnostic.m =================================================================== --- /issm/trunk/src/m/solutions/ice/diagnostic.m (revision 326) +++ /issm/trunk/src/m/solutions/ice/diagnostic.m (revision 327) @@ -54,5 +54,4 @@ %First, build elements,grids,loads, etc ... for horizontal, base vertical and vertical model - fem.m_dbv=CreateFemModel(md,'diagnostic_basevert'); fem.m_dv=CreateFemModel(md,'diagnostic_vert'); Index: /issm/trunk/src/m/solutions/ice/diagnostic3d.m =================================================================== --- /issm/trunk/src/m/solutions/ice/diagnostic3d.m (revision 326) +++ /issm/trunk/src/m/solutions/ice/diagnostic3d.m (revision 327) @@ -60,12 +60,7 @@ inputs.velocity_average=velocity_average; -%Compute basal vertical velocity, in 3d -disp(sprintf('\n%s',['computing basal vertical velocity...'])); -u_g_basevert=diagnostic_core_linear(fem.m_dbv,'diagnostic_basevert',inputs); - %update dirichlet boundary conditions with these new base vertical velocities m_dv=fem.m_dv; gridset=m_dv.gridset; -m_dv.y_g=u_g_basevert; m_dv.ys=Reducevector_g(m_dv.y_g); %compute vertical velocities Index: /issm/trunk/src/m/solutions/ice/diagnostic_core_linear.m =================================================================== --- /issm/trunk/src/m/solutions/ice/diagnostic_core_linear.m (revision 326) +++ /issm/trunk/src/m/solutions/ice/diagnostic_core_linear.m (revision 327) @@ -44,4 +44,9 @@ [K_gg,p_g]=PenaltySystemMatrices(grids,loads,materials,kflag, pflag, sparsity,inputs,analysis_type,K_gg,p_g); +if strcmpi(analysis_type,'diagnostic_hutter'), + K_gg + p_g + %error +end if cluster, K_gg=distributed(gplus(K_gg),'convert');