Changeset 5433

Timestamp:

08/19/10 17:00:42 (15 years ago)

Author:

seroussi

Message:

stokes now independent

Location:

issm/trunk/src

Files:

• c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (modified) (1 diff)
• c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp (modified) (1 diff)
• c/modules/ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp (modified) (3 diffs)
• c/modules/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp (modified) (1 diff)
• c/objects/Elements/Penta.cpp (modified) (9 diffs)
• c/objects/Node.cpp (modified) (2 diffs)
• c/solutions/diagnostic_core.cpp (modified) (1 diff)
• m/classes/public/collapse.m (modified) (1 diff)
• m/classes/public/extrude.m (modified) (3 diffs)
• m/classes/public/ismodelselfconsistent.m (modified) (3 diffs)
• m/classes/public/mechanicalproperties.m (modified) (1 diff)
• m/classes/public/modelextract.m (modified) (1 diff)
• m/classes/public/plot/plot_elementstype.m (modified) (7 diffs)
• m/classes/public/setelementstype.m (modified) (8 diffs)
• m/solutions/SolutionConfiguration.m (modified) (1 diff)
• m/solutions/diagnostic_core.m (modified) (1 diff)

issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp

@@ -46 +46 @@
 
                         /*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]){

issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp

@@ -56 +56 @@
                 
                 /*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++;

issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp

@@ -32 +32 @@
         /*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 +47 @@
 
                                 /*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 +77 @@
 
         /*Free data: */
+        xfree((void**)&iomodel->spcvelocity);
         xfree((void**)&iomodel->gridonstokes);
+        xfree((void**)&iomodel->surface);
+        xfree((void**)&iomodel->z);
 
         cleanup_and_return:

issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp

@@ -56 +56 @@
         
                 /*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;
                 }

issm/trunk/src/c/objects/Elements/Penta.cpp

@@ -2158 +2158 @@
                         }
                         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 +2178 @@
         //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 +2204 @@
         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 +2876 @@
 
         /*Material properties: */
-        double         gravity,rho_ice,rho_water;
+        double    gravity,rho_ice,rho_water;
 
         /*Collapsed formulation: */
@@ -2925 +2933 @@
         double stokesreconditioning;
         int analysis_type;
+        int approximation;
 
         /*inputs: */
@@ -2945 +2954 @@
         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 +4127 @@
         bool shelf;
         bool isstokes;
+        int  approximation;
         Input* vx_input=NULL;
         Input* vy_input=NULL;
@@ -4127 +4138 @@
         inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
         inputs->GetParameterValue(&isstokes,IsStokesEnum);
+        inputs->GetParameterValue(&approximation,ApproximationEnum);
 
         /*retrieve some parameters: */
@@ -4132 +4144 @@
 
         /*If on water or not Stokes, skip load: */
-        if(onwater || !isstokes) return;
+        if(onwater || !isstokes || approximation!=StokesApproximationEnum) return;
 
         /*If on water, skip load: */

issm/trunk/src/c/objects/Node.cpp

@@ -71 +71 @@
 
         /*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 +107 @@
         }
 
-        /*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){

issm/trunk/src/c/solutions/diagnostic_core.cpp

@@ -71 +71 @@
         
         
-        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);
         }
 

issm/trunk/src/m/classes/public/collapse.m

@@ -52 +52 @@
         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);

issm/trunk/src/m/classes/public/extrude.m

@@ -176 +176 @@
 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 +185 @@
         %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 +192 @@
 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
 

issm/trunk/src/m/classes/public/ismodelselfconsistent.m

@@ -56 +56 @@
 %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 +69 @@
 %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 +201 @@
 
         %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

issm/trunk/src/m/classes/public/mechanicalproperties.m

@@ -20 +20 @@
         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

issm/trunk/src/m/classes/public/modelextract.m

@@ -139 +139 @@
                 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);

issm/trunk/src/m/classes/public/plot/plot_elementstype.m

@@ -9 +9 @@
 %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 +19 @@
 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 +46 @@
 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 +54 @@
         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 +62 @@
         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 +70 @@
         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 +79 @@
         %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);

issm/trunk/src/m/classes/public/setelementstype.m

@@ -47 +47 @@
 
 %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 +66 @@
 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 +86 @@
 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 +92 @@
 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 +98 @@
 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 +145 @@
 
         %Create MacaAyealPattynApproximation where needed
-        md.elements_type(find(macayealpattynflag),1)=MacAyealPattynApproximationEnum();
+        md.elements_type(find(macayealpattynflag))=MacAyealPattynApproximationEnum();
 
         %Now recreate gridonmacayeal and gridonmacayealpattyn
@@ -148 +153 @@
 
 %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 +180 @@
 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

issm/trunk/src/m/solutions/SolutionConfiguration.m

@@ -12 +12 @@
         case DiagnosticSolutionEnum,
                 numanalyses=6;
+                %analyses=[DiagnosticHorizAnalysisEnum];
                 analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticStokesAnalysisEnum;DiagnosticHutterAnalysisEnum;SurfaceSlopeAnalysisEnum;BedSlopeAnalysisEnum];
 

issm/trunk/src/m/solutions/diagnostic_core.m

@@ -55 +55 @@
         
 
-        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