Changeset 19291

Timestamp:

04/21/15 10:23:49 (10 years ago)

Author:

Mathieu Morlighem

Message:

BUG: fixed spc transient in SIA (if no coupling)

File:

• issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp (modified) (1 diff)

issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp

@@ -9 +9 @@
 void StressbalanceSIAAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
 
+        /*Intermediaries*/
+        bool       isSIA,isSSA,isL1L2,isHO,isFS,iscoupling;
+
         /*Fetch parameters: */
-        bool isSIA;
         iomodel->Constant(&isSIA,FlowequationIsSIAEnum);
+        iomodel->Constant(&isSSA,FlowequationIsSSAEnum);
+        iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum);
+        iomodel->Constant(&isHO,FlowequationIsHOEnum);
+        iomodel->Constant(&isFS,FlowequationIsFSEnum);
 
         /*Now, is the flag isSIA on? otherwise, do nothing: */
         if (!isSIA) return;
 
-        /*Fetch data: */
-        iomodel->FetchData(3,StressbalanceSpcvxEnum,StressbalanceSpcvyEnum,FlowequationVertexEquationEnum);
-
-        /*Initialize conunter*/
-        int count=0;
-
-        /*vx and vy are spc'd if we are not on nodeonSIA: */
-        for(int i=0;i<iomodel->numberofvertices;i++){
-                /*keep only this partition's nodes:*/
-                if((iomodel->my_vertices[i])){
-                        if (reCast<int,IssmDouble>(iomodel->Data(FlowequationVertexEquationEnum)[i])!=SIAApproximationEnum){
-
-                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,0,StressbalanceSIAAnalysisEnum));
-                                count++;
-
-                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceSIAAnalysisEnum));
-                                count++;
-                        }
-                        else{
-                                if (!xIsNan<IssmDouble>(iomodel->Data(StressbalanceSpcvxEnum)[i])){
-                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,iomodel->Data(StressbalanceSpcvxEnum)[i],StressbalanceSIAAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
+        /*Do we have coupling*/
+        if((isSIA?1.:0.) + (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
+         iscoupling = true;
+        else
+         iscoupling = false;
+
+        /*If no coupling, call Regular IoModelToConstraintsx, else, OLD stuff, keep for now*/
+        if(!iscoupling){
+                IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceSIAAnalysisEnum,P1Enum,0);
+                IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvyEnum,StressbalanceSIAAnalysisEnum,P1Enum,1);
+        }
+        else{
+                /*Fetch data: */
+                iomodel->FetchData(3,StressbalanceSpcvxEnum,StressbalanceSpcvyEnum,FlowequationVertexEquationEnum);
+
+                /*Initialize conunter*/
+                int count=0;
+
+                /*vx and vy are spc'd if we are not on nodeonSIA: */
+                for(int i=0;i<iomodel->numberofvertices;i++){
+                        /*keep only this partition's nodes:*/
+                        if((iomodel->my_vertices[i])){
+                                if (reCast<int,IssmDouble>(iomodel->Data(FlowequationVertexEquationEnum)[i])!=SIAApproximationEnum){
+
+                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,0,StressbalanceSIAAnalysisEnum));
+                                        count++;
+
+                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceSIAAnalysisEnum));
                                         count++;
                                 }
-
-                                if (!xIsNan<IssmDouble>(iomodel->Data(StressbalanceSpcvyEnum)[i])){
-                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,iomodel->Data(StressbalanceSpcvyEnum)[i],StressbalanceSIAAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy
-                                        count++;
+                                else{
+                                        if (!xIsNan<IssmDouble>(iomodel->Data(StressbalanceSpcvxEnum)[i])){
+                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,iomodel->Data(StressbalanceSpcvxEnum)[i],StressbalanceSIAAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
+                                                count++;
+                                        }
+
+                                        if (!xIsNan<IssmDouble>(iomodel->Data(StressbalanceSpcvyEnum)[i])){
+                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,iomodel->Data(StressbalanceSpcvyEnum)[i],StressbalanceSIAAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy
+                                                count++;
+                                        }
                                 }
                         }