Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (revision 5432) +++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (revision 5433) @@ -46,5 +46,5 @@ /*Start with adding spcs of coupling: zero at the border macayeal/pattyn for the appropriate dofs*/ - if ((int)iomodel->vertices_type[2*i+0]==MacAyealPattynApproximationEnum){ + if ((int)iomodel->vertices_type[i]==MacAyealPattynApproximationEnum){ /*If grionmacayeal, spc pattyn dofs: 3 & 4*/ if ((int)iomodel->gridonpattyn[i]){ Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp (revision 5432) +++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp (revision 5433) @@ -56,16 +56,16 @@ /*Do not create ice front if Hutter or Stokes elements*/ - if ((int)*(iomodel->elements_type+2*element+0)==(HutterApproximationEnum || StokesApproximationEnum)) continue; + if ((int)*(iomodel->elements_type+element)==(HutterApproximationEnum || StokesApproximationEnum)) continue; /*Create and add load: */ - if ((int)*(iomodel->elements_type+2*element+0)==(MacAyealApproximationEnum)){ + if ((int)*(iomodel->elements_type+element)==(MacAyealApproximationEnum)){ loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,SegmentIcefrontEnum,DiagnosticHorizAnalysisEnum)); count++; } - else if ((int)*(iomodel->elements_type+2*element+0)==(PattynApproximationEnum)){ + else if ((int)*(iomodel->elements_type+element)==(PattynApproximationEnum)){ loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,QuadIceFrontEnum,DiagnosticHorizAnalysisEnum)); count++; } - else if ((int)*(iomodel->elements_type+2*element+0)==(MacAyealPattynApproximationEnum)){ + else if ((int)*(iomodel->elements_type+element)==(MacAyealPattynApproximationEnum)){ loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,SegmentIcefrontEnum,DiagnosticHorizAnalysisEnum)); count++; Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp (revision 5432) +++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp (revision 5433) @@ -32,4 +32,7 @@ /*Fetch data: */ IoModelFetchData(&iomodel->gridonstokes,NULL,NULL,iomodel_handle,"gridonstokes"); + IoModelFetchData(&iomodel->spcvelocity,NULL,NULL,iomodel_handle,"spcvelocity"); + IoModelFetchData(&iomodel->surface,NULL,NULL,iomodel_handle,"surface"); + IoModelFetchData(&iomodel->z,NULL,NULL,iomodel_handle,"z"); /*Initialize counter*/ @@ -44,12 +47,29 @@ /*This grid will see its vx,vy and vz dofs spc'd to pattyn velocities: */ - constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,DiagnosticStokesAnalysisEnum)); + constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,*(iomodel->spcvelocity+6*i+3)/iomodel->yts,DiagnosticStokesAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx count++; - constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,DiagnosticStokesAnalysisEnum)); + constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,*(iomodel->spcvelocity+6*i+4)/iomodel->yts,DiagnosticStokesAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vy count++; - constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,0,DiagnosticStokesAnalysisEnum)); + constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,*(iomodel->spcvelocity+6*i+5)/iomodel->yts,DiagnosticStokesAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vz count++; + + constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,iomodel->g*iomodel->rho_ice*(iomodel->surface[i]-iomodel->z[i])/iomodel->stokesreconditioning,DiagnosticStokesAnalysisEnum)); //add lithostatic spc for pressure + count++; + } + else{ //this is a regular node so use spcs if needed + if ((int)iomodel->spcvelocity[6*i+0]){ + constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,*(iomodel->spcvelocity+6*i+3)/iomodel->yts,DiagnosticStokesAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + } + if ((int)iomodel->spcvelocity[6*i+1]){ + constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,*(iomodel->spcvelocity+6*i+4)/iomodel->yts,DiagnosticStokesAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy + count++; + } + if ((int)iomodel->spcvelocity[6*i+2]){ + constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,*(iomodel->spcvelocity+6*i+5)/iomodel->yts,DiagnosticStokesAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vz + count++; + } } } //if((my_vertices[i]==1)) @@ -57,5 +77,8 @@ /*Free data: */ + xfree((void**)&iomodel->spcvelocity); xfree((void**)&iomodel->gridonstokes); + xfree((void**)&iomodel->surface); + xfree((void**)&iomodel->z); cleanup_and_return: Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp (revision 5432) +++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp (revision 5433) @@ -56,5 +56,5 @@ /*Do not create ice front if it is not a Stokes element*/ - if ((int)*(iomodel->elements_type+2*element+1)!=StokesApproximationEnum){ + if ((int)*(iomodel->elements_type+element)!=StokesApproximationEnum){ continue; } Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5432) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5433) @@ -2158,9 +2158,13 @@ } if(!iomodel->pressure){ - for(i=0;i<6;i++)nodeinputs[i]=0; + for(i=0;i<6;i++)nodeinputs[i]=i; if(iomodel->qmu_analysis){ this->inputs->AddInput(new PentaVertexInput(PressureEnum,nodeinputs)); this->inputs->AddInput(new PentaVertexInput(QmuPressureEnum,nodeinputs)); } + if(iomodel->isstokes){ + this->inputs->AddInput(new PentaVertexInput(PressureEnum,nodeinputs)); + this->inputs->AddInput(new PentaVertexInput(PressureOldEnum,nodeinputs)); + } } break; @@ -2174,20 +2178,24 @@ //Need to know the type of approximation for this element if(iomodel->elements_type){ - if (*(iomodel->elements_type+2*index+0)==MacAyealApproximationEnum){ + if (*(iomodel->elements_type+index)==MacAyealApproximationEnum){ this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealApproximationEnum)); } + else if (*(iomodel->elements_type+index)==PattynApproximationEnum){ + this->inputs->AddInput(new IntInput(ApproximationEnum,PattynApproximationEnum)); + } + else if (*(iomodel->elements_type+index)==MacAyealPattynApproximationEnum){ + this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealPattynApproximationEnum)); + } + else if (*(iomodel->elements_type+index)==HutterApproximationEnum){ + this->inputs->AddInput(new IntInput(ApproximationEnum,HutterApproximationEnum)); + } + else if (*(iomodel->elements_type+index)==StokesApproximationEnum){ + this->inputs->AddInput(new IntInput(ApproximationEnum,StokesApproximationEnum)); + } + else if (*(iomodel->elements_type+index)==NoneApproximationEnum){ + this->inputs->AddInput(new IntInput(ApproximationEnum,NoneApproximationEnum)); + } else{ - if (*(iomodel->elements_type+2*index+0)==PattynApproximationEnum){ - this->inputs->AddInput(new IntInput(ApproximationEnum,PattynApproximationEnum)); - } - else if (*(iomodel->elements_type+2*index+0)==MacAyealPattynApproximationEnum){ - this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealPattynApproximationEnum)); - } - else if (*(iomodel->elements_type+2*index+0)==HutterApproximationEnum){ - this->inputs->AddInput(new IntInput(ApproximationEnum,HutterApproximationEnum)); - } - else{ - ISSMERROR("Approximation type %s not supported yet",EnumToString((int)*(iomodel->elements_type+2*index+0))); - } + ISSMERROR("Approximation type %s not supported yet",EnumToString((int)*(iomodel->elements_type+index))); } } @@ -2196,5 +2204,5 @@ if(analysis_type==DiagnosticStokesAnalysisEnum){ ISSMASSERT(iomodel->elements_type); - if(iomodel->elements_type[2*index+1]==StokesApproximationEnum) + if(iomodel->elements_type[index]==StokesApproximationEnum) this->inputs->AddInput(new BoolInput(IsStokesEnum,true)); else @@ -2868,5 +2876,5 @@ /*Material properties: */ - double gravity,rho_ice,rho_water; + double gravity,rho_ice,rho_water; /*Collapsed formulation: */ @@ -2925,4 +2933,5 @@ double stokesreconditioning; int analysis_type; + int approximation; /*inputs: */ @@ -2945,7 +2954,8 @@ inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); inputs->GetParameterValue(&isstokes,IsStokesEnum); + inputs->GetParameterValue(&approximation,ApproximationEnum); /*If on water or not Stokes, skip stiffness: */ - if(onwater || !isstokes) return; + if(onwater || !isstokes || approximation!=StokesApproximationEnum) return; /*recovre material parameters: */ @@ -4117,4 +4127,5 @@ bool shelf; bool isstokes; + int approximation; Input* vx_input=NULL; Input* vy_input=NULL; @@ -4127,4 +4138,5 @@ inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); inputs->GetParameterValue(&isstokes,IsStokesEnum); + inputs->GetParameterValue(&approximation,ApproximationEnum); /*retrieve some parameters: */ @@ -4132,5 +4144,5 @@ /*If on water or not Stokes, skip load: */ - if(onwater || !isstokes) return; + if(onwater || !isstokes || approximation!=StokesApproximationEnum) return; /*If on water, skip load: */ Index: /issm/trunk/src/c/objects/Node.cpp =================================================================== --- /issm/trunk/src/c/objects/Node.cpp (revision 5432) +++ /issm/trunk/src/c/objects/Node.cpp (revision 5433) @@ -71,5 +71,5 @@ /*indexing:*/ - DistributeNumDofs(&this->indexing,analysis_type,iomodel->vertices_type+2*io_index); //number of dofs per node + DistributeNumDofs(&this->indexing,analysis_type,iomodel->vertices_type+io_index); //number of dofs per node numdofs=this->indexing.numberofdofs; @@ -107,21 +107,4 @@ } - /*Diagnostic Stokes*/ - if (analysis_type==DiagnosticStokesAnalysisEnum){ - /*On a 3d mesh, in stokes formualtions, only stokes grids are free, the others are frozen: */ - if (!iomodel->borderstokes) ISSMERROR("iomodel->borderstokes is NULL"); - if (iomodel->borderstokes[io_index]){ - //freeze everything except pressure - this->FreezeDof(1); - this->FreezeDof(2); - this->FreezeDof(3); - } - else if (iomodel->gridonstokes[io_index]==0){ - for(k=1;k<=numdofs;k++){ - this->FreezeDof(k); - } - } - } - /*Diagnostic Hutter*/ if (analysis_type==DiagnosticHutterAnalysisEnum){ Index: /issm/trunk/src/c/solutions/diagnostic_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/diagnostic_core.cpp (revision 5432) +++ /issm/trunk/src/c/solutions/diagnostic_core.cpp (revision 5433) @@ -71,19 +71,19 @@ - if (dim==3){ + if (dim==3 & ~isstokes){ if(verbose)_printf_("%s\n"," computing vertical velocities..."); femmodel->SetCurrentConfiguration(DiagnosticVertAnalysisEnum); solver_linear(femmodel); + } - if (isstokes){ + if (isstokes){ - if(verbose)_printf_("%s\n"," update boundary conditions for stokes using velocities previously computed..."); - ResetBoundaryConditions(femmodel,DiagnosticStokesAnalysisEnum); + // if(verbose)_printf_("%s\n"," update boundary conditions for stokes using velocities previously computed..."); + // ResetBoundaryConditions(femmodel,DiagnosticStokesAnalysisEnum); - if(verbose)_printf_("%s\n"," computing stokes velocities and pressure ..."); - femmodel->SetCurrentConfiguration(DiagnosticStokesAnalysisEnum); - solver_diagnostic_nonlinear(femmodel,conserve_loads); - } + if(verbose)_printf_("%s\n"," computing stokes velocities and pressure ..."); + femmodel->SetCurrentConfiguration(DiagnosticStokesAnalysisEnum); + solver_diagnostic_nonlinear(femmodel,conserve_loads); } Index: /issm/trunk/src/m/classes/public/collapse.m =================================================================== --- /issm/trunk/src/m/classes/public/collapse.m (revision 5432) +++ /issm/trunk/src/m/classes/public/collapse.m (revision 5433) @@ -52,5 +52,5 @@ md.elements_type=md.elements_type2d; elseif ~isnan(md.elements_type) - md.elements_type=project2d(md,md.elements_type,1); + md.elements_type=project2d(md,md.elements_type); end md.gridonhutter=project2d(md,md.gridonhutter,1); Index: /issm/trunk/src/m/classes/public/extrude.m =================================================================== --- /issm/trunk/src/m/classes/public/extrude.m (revision 5432) +++ /issm/trunk/src/m/classes/public/extrude.m (revision 5433) @@ -176,7 +176,6 @@ if ~isnan(md.elements_type) oldelements_type=md.elements_type2d; - md.elements_type=zeros(number_el3d,2); - md.elements_type(:,1)=project3d(md,oldelements_type(:,1),'element'); - md.elements_type(:,2)=project3d(md,oldelements_type(:,2),'element'); + md.elements_type=zeros(number_el3d,1); + md.elements_type=project3d(md,oldelements_type,'element'); md.gridonhutter=project3d(md,md.gridonhutter,'node'); md.gridonmacayeal=project3d(md,md.gridonmacayeal,'node'); @@ -186,5 +185,5 @@ %dead grids md.deadgrids=ones(md.numberofgrids,1); - md.deadgrids(md.elements(md.elements_type(:,1)~=MacAyealApproximationEnum,:))=0;%non macayeal grids are not dead + md.deadgrids(md.elements(md.elements_type~=MacAyealApproximationEnum,:))=0;%non macayeal grids are not dead md.deadgrids(find(md.gridonbed))=0;%grids from elements on bed are not dead end @@ -193,7 +192,6 @@ if ~isnan(md.vertices_type) oldvertices_type=md.vertices_type2d; - md.vertices_type=zeros(number_grids3d,2); - md.vertices_type(:,1)=project3d(md,oldvertices_type(:,1),'node'); - md.vertices_type(:,2)=project3d(md,oldvertices_type(:,2),'node'); + md.vertices_type=zeros(number_grids3d,1); + md.vertices_type=project3d(md,oldvertices_type,'node'); end Index: /issm/trunk/src/m/classes/public/ismodelselfconsistent.m =================================================================== --- /issm/trunk/src/m/classes/public/ismodelselfconsistent.m (revision 5432) +++ /issm/trunk/src/m/classes/public/ismodelselfconsistent.m (revision 5433) @@ -56,11 +56,9 @@ %Check the size of elements_type fields={'elements_type'}; -checksize(md,fields,[md.numberofelements 2]); -%Check the values of elements_type(1) -checkvalues(md,{'elements_type(:,1)'},[MacAyealApproximationEnum() HutterApproximationEnum() PattynApproximationEnum() MacAyealPattynApproximationEnum()]); -%Check the values of elements_type(2) -checkvalues(md,{'elements_type(:,2)'},[StokesApproximationEnum() NoneApproximationEnum()]); +checksize(md,fields,[md.numberofelements 1]); +%Check the values of elements_type +checkvalues(md,{'elements_type'},[MacAyealApproximationEnum() HutterApproximationEnum() PattynApproximationEnum() MacAyealPattynApproximationEnum() StokesApproximationEnum() NoneApproximationEnum()]); if (md.dim==2), - checkvalues(md,{'elements_type(:,1)'},[MacAyealApproximationEnum() HutterApproximationEnum()]); + checkvalues(md,{'elements_type'},[MacAyealApproximationEnum() HutterApproximationEnum()]); end if (md.ismacayealpattyn==0 && md.ishutter==0 && md.isstokes==0), @@ -71,11 +69,9 @@ %Check the size of verticess_type fields={'vertices_type'}; -checksize(md,fields,[md.numberofgrids 2]); -%Check the values of elements_type(1) -checkvalues(md,{'elements_type(:,1)'},[MacAyealApproximationEnum() HutterApproximationEnum() PattynApproximationEnum() MacAyealPattynApproximationEnum()]); -%Check the values of elements_type(2) -checkvalues(md,{'elements_type(:,2)'},[StokesApproximationEnum() NoneApproximationEnum()]); +checksize(md,fields,[md.numberofgrids 1]); +%Check the values of vertices_type +checkvalues(md,{'vertices_type'},[MacAyealApproximationEnum() HutterApproximationEnum() PattynApproximationEnum() MacAyealPattynApproximationEnum() StokesApproximationEnum() NoneApproximationEnum()]); if (md.dim==2), - checkvalues(md,{'elements_type(:,1)'},[MacAyealApproximationEnum() HutterApproximationEnum()]); + checkvalues(md,{'vertices_type'},[MacAyealApproximationEnum() HutterApproximationEnum()]); end if (md.ismacayealpattyn==0 && md.ishutter==0 && md.isstokes==0), @@ -205,5 +201,5 @@ %HUTTER ON ICESHELF WARNING - if any(md.elements_type(:,1)==HutterApproximationEnum & md.elementoniceshelf), + if any(md.elements_type==HutterApproximationEnum & md.elementoniceshelf), disp(sprintf('\n !!! Warning: Hutter''s model is not consistent on ice shelves !!!\n')); end Index: /issm/trunk/src/m/classes/public/mechanicalproperties.m =================================================================== --- /issm/trunk/src/m/classes/public/mechanicalproperties.m (revision 5432) +++ /issm/trunk/src/m/classes/public/mechanicalproperties.m (revision 5433) @@ -20,5 +20,5 @@ error('only 2d model supported yet'); end -if any(md.elements_type(:,1)~=MacAyealApproximationEnum), +if any(md.elements_type~=MacAyealApproximationEnum), disp('Warning: the model has some non macayeal elements. These will be treated like MacAyeal''s elements'); end Index: /issm/trunk/src/m/classes/public/modelextract.m =================================================================== --- /issm/trunk/src/m/classes/public/modelextract.m (revision 5432) +++ /issm/trunk/src/m/classes/public/modelextract.m (revision 5433) @@ -139,5 +139,5 @@ md2.elements2d(:,3)=Pgrid(md2.elements2d(:,3)); - if ~isnan(md2.elements_type2d), md2.elements_type2d=md1.elements_type2d(pos_elem_2d,:); end; + if ~isnan(md2.elements_type2d), md2.elements_type2d=md1.elements_type2d(pos_elem_2d); end; md2.x2d=md1.x(pos_grid_2d); md2.y2d=md1.y(pos_grid_2d); Index: /issm/trunk/src/m/classes/public/plot/plot_elementstype.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plot_elementstype.m (revision 5432) +++ /issm/trunk/src/m/classes/public/plot/plot_elementstype.m (revision 5433) @@ -9,7 +9,5 @@ %process data and model [x y z elements is2d]=processmesh(md,[],options); -[data1 datatype]=processdata(md,md.elements_type(:,1),options); -[data2 datatype]=processdata(md,md.elements_type(:,2),options); -data=[data1 data2]; +[data datatype]=processdata(md,md.elements_type,options); %edgecolor? @@ -21,22 +19,22 @@ if is2d %Hutter elements - posH=find(data(:,1)==HutterApproximationEnum); + posH=find(data==HutterApproximationEnum); A=elements(posH,1); B=elements(posH,2); C=elements(posH,3); p1=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',HutterApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor); %MacAyeal element - posH=find(data(:,1)==MacAyealApproximationEnum); + posH=find(data==MacAyealApproximationEnum); A=elements(posH,1); B=elements(posH,2); C=elements(posH,3); p2=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',MacAyealApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor); %Pattyn element - posH=find(data(:,1)==PattynApproximationEnum); + posH=find(data==PattynApproximationEnum); A=elements(posH,1); B=elements(posH,2); C=elements(posH,3); p3=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',PattynApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor); %MacAyealPattyn element - posH=find(data(:,1)==MacAyealPattynApproximationEnum); + posH=find(data==MacAyealPattynApproximationEnum); A=elements(posH,1); B=elements(posH,2); C=elements(posH,3); p5=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',MacAyealPattynApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor); %Stokes elements alpha=0.35; - posS=find(data(:,2)==StokesApproximationEnum); + posS=find(data==StokesApproximationEnum); if ~isempty(posS) A=elements(posS,1); B=elements(posS,2); C=elements(posS,3); @@ -48,5 +46,5 @@ else %Hutter elements - posH=find(data(:,1)==HutterApproximationEnum); + posH=find(data==HutterApproximationEnum); A=elements(posH,1); B=elements(posH,2); C=elements(posH,3); D=elements(posH,4); E=elements(posH,5); F=elements(posH,6); p1=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData', HutterApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor); @@ -56,5 +54,5 @@ patch( 'Faces', [C A D F ], 'Vertices', [x y z],'CData', HutterApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor); %MacAyeal elements - posM=find(data(:,1)==MacAyealApproximationEnum); + posM=find(data==MacAyealApproximationEnum); A=elements(posM,1); B=elements(posM,2); C=elements(posM,3); D=elements(posM,4); E=elements(posM,5); F=elements(posM,6); p2=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData', MacAyealApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor); @@ -64,5 +62,5 @@ patch( 'Faces', [C A D F ], 'Vertices', [x y z],'CData', MacAyealApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor); %Pattyn elements - posP=find(data(:,1)==PattynApproximationEnum); + posP=find(data==PattynApproximationEnum); A=elements(posP,1); B=elements(posP,2); C=elements(posP,3); D=elements(posP,4); E=elements(posP,5); F=elements(posP,6); p3=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData', PattynApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor); @@ -72,5 +70,5 @@ patch( 'Faces', [C A D F ], 'Vertices', [x y z],'CData', PattynApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor); %MacAyealPattyn elements - posP=find(data(:,1)==MacAyealPattynApproximationEnum); + posP=find(data==MacAyealPattynApproximationEnum); A=elements(posP,1); B=elements(posP,2); C=elements(posP,3); D=elements(posP,4); E=elements(posP,5); F=elements(posP,6); p5=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData', MacAyealPattynApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor); @@ -81,5 +79,5 @@ %Stokes elements alpha=0.35; - posS=find(data(:,2)==StokesApproximationEnum); + posS=find(data==StokesApproximationEnum); if ~isempty(posS) A=elements(posS,1); B=elements(posS,2); C=elements(posS,3); D=elements(posS,4); E=elements(posS,5); F=elements(posS,6); Index: /issm/trunk/src/m/classes/public/setelementstype.m =================================================================== --- /issm/trunk/src/m/classes/public/setelementstype.m (revision 5432) +++ /issm/trunk/src/m/classes/public/setelementstype.m (revision 5433) @@ -47,5 +47,5 @@ %check that each element has at least one flag -if any(hutterflag+ macayealflag+pattynflag==0), +if any(hutterflag+ macayealflag+pattynflag+stokesflag==0), error('setelementstype error message: elements type not assigned, must be specified') end @@ -66,19 +66,19 @@ end -%Modify stokesflag to get rid of elements contrained everywhere -nonstokesflag=find(~stokesflag); %non stokes elements -gridonstokes=ones(md.numberofgrids,1); -gridonstokes(md.elements(nonstokesflag,:))=0; %non stokes grids -pos=find(sum(md.spcvelocity(:,1:2),2)==2); %find all the grids on the boundary of the domain without icefront -gridonstokes(pos)=0; %we have to constrain all the boundary of the domain without icefront +%Stokes can only be used alone for now: +if any(stokesflag) &any(hutterflag+pattynflag+macayealflag), + error('setelementstype error message: stokes cannot be used with any other model for now, put stokes everywhere') +end -%Rule out elements where all grids are on the border of the stokes grids (because those elements are completely -%constrained -pos=find(stokesflag); -constrained_list=pos(find(sum(gridonstokes(md.elements(pos,:)),2)==0)); -stokesflag(constrained_list)=0; - +if any(stokesflag), + %Modify stokesflag to get rid of elements contrained everywhere + fullspcnodes=double(sum(md.spcvelocity(:,1:3),2)==3); %find all the grids on the boundary of the domain without icefront + fullspcelems=double(sum(fullspcnodes(md.elements),2)==6); %find all the grids on the boundary of the domain without icefront + stokesflag(find(fullspcelems))=0; + gridonstokes=zeros(md.numberofgrids,1); +end + %add in model who is who -md.elements_type=zeros(md.numberofelements,2); +md.elements_type=zeros(md.numberofelements,1); %1: Hutter elements @@ -86,5 +86,5 @@ gridonhutter(md.elements(find(hutterflag),:))=1; md.gridonhutter=gridonhutter; -md.elements_type(find(hutterflag),1)=HutterApproximationEnum(); +md.elements_type(find(hutterflag))=HutterApproximationEnum(); %2: MacAyeal elements @@ -92,5 +92,5 @@ gridonmacayeal(md.elements(find(macayealflag),:))=1; md.gridonmacayeal=gridonmacayeal; -md.elements_type(find(macayealflag),1)=MacAyealApproximationEnum(); +md.elements_type(find(macayealflag))=MacAyealApproximationEnum(); %3: Pattyn elements @@ -98,12 +98,17 @@ gridonpattyn(md.elements(find(pattynflag),:))=1; md.gridonpattyn=gridonpattyn; -md.elements_type(find(pattynflag),1)=PattynApproximationEnum(); +md.elements_type(find(pattynflag))=PattynApproximationEnum(); %4: Stokes elements +gridonstokes=zeros(md.numberofgrids,1); +gridonstokes(md.elements(find(stokesflag),:))=1; md.gridonstokes=gridonstokes; -md.elements_type(find(stokesflag),2)=StokesApproximationEnum(); +md.elements_type(find(stokesflag))=StokesApproximationEnum(); -%5: None elements (non Stokes) -md.elements_type(find(~stokesflag),2)=NoneApproximationEnum(); +%Complete with NoneApproximation if there is no stokes +if any(stokesflag), + %5: None elements (non Stokes) + md.elements_type(find(~stokesflag))=NoneApproximationEnum(); +end %Now take care of the coupling between MacAyeal and Pattyn @@ -140,5 +145,5 @@ %Create MacaAyealPattynApproximation where needed - md.elements_type(find(macayealpattynflag),1)=MacAyealPattynApproximationEnum(); + md.elements_type(find(macayealpattynflag))=MacAyealPattynApproximationEnum(); %Now recreate gridonmacayeal and gridonmacayealpattyn @@ -148,21 +153,23 @@ %Create vertices_type -md.vertices_type=zeros(md.numberofgrids,2); +md.vertices_type=zeros(md.numberofgrids,1); pos=find(gridonhutter); -md.vertices_type(pos,1)=HutterApproximationEnum(); +md.vertices_type(pos)=HutterApproximationEnum(); pos=find(gridonmacayeal); -md.vertices_type(pos,1)=MacAyealApproximationEnum(); +md.vertices_type(pos)=MacAyealApproximationEnum(); pos=find(gridonpattyn); -md.vertices_type(pos,1)=PattynApproximationEnum(); +md.vertices_type(pos)=PattynApproximationEnum(); pos=find(gridonhutter); -md.vertices_type(pos,1)=HutterApproximationEnum(); +md.vertices_type(pos)=HutterApproximationEnum(); pos=find(gridonpattyn & gridonmacayeal); -md.vertices_type(pos,1)=PattynApproximationEnum(); +md.vertices_type(pos)=PattynApproximationEnum(); pos=find(gridonmacayealpattyn); -md.vertices_type(pos,1)=MacAyealPattynApproximationEnum(); +md.vertices_type(pos)=MacAyealPattynApproximationEnum(); pos=find(gridonstokes); -md.vertices_type(pos,2)=StokesApproximationEnum(); -pos=find(~gridonstokes); -md.vertices_type(pos,2)=NoneApproximationEnum(); +md.vertices_type(pos)=StokesApproximationEnum(); +if any(stokesflag), + pos=find(~gridonstokes); + md.vertices_type(pos)=NoneApproximationEnum(); +end %flag dead grids (strictly in MacAyeal and not on bed -> not used in diagnostic horiz) @@ -173,14 +180,8 @@ md.deadgrids=deadgrids; -%figure out the border stokes grids -stokes_elements=find(md.elements_type(:,2)==StokesApproximationEnum()); %find the elements on the stokes domain -borderflags=zeros(md.numberofgrids,1); -borderflags(md.elements(stokes_elements,:))=1; %find all the grids of the elements on stokes domain, ie stokes grids and borderstokes -md.borderstokes=borderflags-md.gridonstokes; %remove stokes grids from this list - %figure out solution types -md.ishutter=double(any(md.elements_type(:,1)==HutterApproximationEnum)); -md.ismacayealpattyn=double(any(md.elements_type(:,1)==MacAyealApproximationEnum | md.elements_type(:,1)==PattynApproximationEnum)); -md.isstokes=double(any(md.elements_type(:,2)==StokesApproximationEnum)); +md.ishutter=double(any(md.elements_type==HutterApproximationEnum)); +md.ismacayealpattyn=double(any(md.elements_type==MacAyealApproximationEnum | md.elements_type==PattynApproximationEnum)); +md.isstokes=double(any(md.elements_type==StokesApproximationEnum)); end Index: /issm/trunk/src/m/solutions/SolutionConfiguration.m =================================================================== --- /issm/trunk/src/m/solutions/SolutionConfiguration.m (revision 5432) +++ /issm/trunk/src/m/solutions/SolutionConfiguration.m (revision 5433) @@ -12,4 +12,5 @@ case DiagnosticSolutionEnum, numanalyses=6; + %analyses=[DiagnosticHorizAnalysisEnum]; analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticStokesAnalysisEnum;DiagnosticHutterAnalysisEnum;SurfaceSlopeAnalysisEnum;BedSlopeAnalysisEnum]; Index: /issm/trunk/src/m/solutions/diagnostic_core.m =================================================================== --- /issm/trunk/src/m/solutions/diagnostic_core.m (revision 5432) +++ /issm/trunk/src/m/solutions/diagnostic_core.m (revision 5433) @@ -55,19 +55,19 @@ - if dim==3, + if dim==3 & ~isstokes, displaystring(verbose,'\n%s',['computing vertical velocities...']); femmodel=SetCurrentConfiguration(femmodel,DiagnosticVertAnalysisEnum); femmodel=solver_linear(femmodel); + end - if isstokes, + if isstokes, - displaystring(verbose,'\n%s',['update boundary conditions for stokes using velocities previously computed...']); - femmodel=ResetBoundaryConditions(femmodel,DiagnosticStokesAnalysisEnum); +% displaystring(verbose,'\n%s',['update boundary conditions for stokes using velocities previously computed...']); +% femmodel=ResetBoundaryConditions(femmodel,DiagnosticStokesAnalysisEnum); - displaystring(verbose,'\n%s',['computing stokes velocities and pressure...']); - femmodel=SetCurrentConfiguration(femmodel,DiagnosticStokesAnalysisEnum); - femmodel=solver_diagnostic_nonlinear(femmodel,conserve_loads); - end + displaystring(verbose,'\n%s',['computing stokes velocities and pressure...']); + femmodel=SetCurrentConfiguration(femmodel,DiagnosticStokesAnalysisEnum); + femmodel=solver_diagnostic_nonlinear(femmodel,conserve_loads); end