Changeset 15624

Timestamp:

07/25/13 15:07:35 (12 years ago)

Author:

Mathieu Morlighem

Message:

NEW: for SSA and HO, use regular IoModelToConstraintx

File:

• issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (modified) (6 diffs)

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

@@ -6 +6 @@
 #include "../../../classes/classes.h"
 #include "../../../shared/shared.h"
+#include "../../IoModelToConstraintsx/IoModelToConstraintsx.h"
 #include "../ModelProcessorx.h"
 
@@ -12 +13 @@
         /*Intermediary*/
         int        i,j;
-        int        count;
+        int        count,temp,finiteelement;
         IssmDouble g;
         IssmDouble rho_ice;
         IssmDouble FSreconditioning;
-        bool       isSSA,isL1L2,isHO,isFS;
-        int        fe_ssa,fe_ho;
+        bool       isSSA,isL1L2,isHO,isFS,iscoupling;
         bool       spcpresent = false;
         int        Mx,Nx;
@@ -50 +50 @@
         iomodel->Constant(&isHO,FlowequationIsHOEnum);
         iomodel->Constant(&isFS,FlowequationIsFSEnum);
-        iomodel->Constant(&fe_ssa,FlowequationFeSSAEnum);
-        iomodel->Constant(&fe_ho,FlowequationFeHOEnum);
 
         /*Recover pointer: */
@@ -58 +56 @@
         /*Now, is the flag macayaealHO on? otherwise, do nothing: */
         if(!isSSA & !isHO & !isFS & !isL1L2){
+                *pconstraints=constraints;
+                return;
+        }
+
+        /*Do we have coupling*/
+        if( (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, use P1 elements only*/
+        if(!iscoupling && !isFS){
+
+                /*Get finite element type*/
+                if(isSSA){
+                        iomodel->Constant(&temp,FlowequationFeSSAEnum);
+                        switch(temp){
+                                case 0 : finiteelement = P1Enum; break;
+                                case 1 : finiteelement = P2Enum; break;
+                                default: _error_("finite element "<<temp<<" not supported");
+                        }
+                }
+                else if(isL1L2){
+                        finiteelement = P1Enum;
+                }
+                else if(isHO){
+                        iomodel->Constant(&temp,FlowequationFeHOEnum);
+                        switch(temp){
+                                case 0 : finiteelement = P1Enum; break;
+                                case 1 : finiteelement = P2Enum; break;
+                                default: _error_("finite element "<<temp<<" not supported");
+                        }
+                }
+                else if(isFS){
+                        finiteelement = P1Enum;
+                }
+                IoModelToConstraintsx(constraints,iomodel,DiagnosticSpcvxEnum,DiagnosticHorizAnalysisEnum,finiteelement,1);
+                IoModelToConstraintsx(constraints,iomodel,DiagnosticSpcvyEnum,DiagnosticHorizAnalysisEnum,finiteelement,2);
+
                 *pconstraints=constraints;
                 return;
@@ -99 +136 @@
 
                         /*Start with adding spcs of coupling: zero at the border SSA/HO for the appropriate dofs*/
-                        if (reCast<int,IssmDouble>(vertices_type[i]==SSAHOApproximationEnum)){
+                        if(reCast<int,IssmDouble>(vertices_type[i]==SSAHOApproximationEnum)){
                                 /*If grionSSA, spc HO dofs: 3 & 4*/
                                         if (reCast<int,IssmDouble>(nodeonHO[i])){
@@ -309 +346 @@
                         }
                 }
-        }
-
-        /*SPC Quadratic elements*/
-        if((isSSA && fe_ssa==1) || (isHO && fe_ho==1)){
-
-                int   v1,v2;
-                bool *my_edges = NULL;
-
-                if(Mx!=iomodel->numberofvertices) _error_("transient spc not supported yet");
-                EdgesPartitioning(&my_edges,iomodel);
-
-                for(i=0;i<iomodel->numberofedges;i++){
-
-                        if(my_edges[i]){
-
-                                v1 = iomodel->edges[2*i+0]-1;
-                                v2 = iomodel->edges[2*i+1]-1;
-
-                                if(!xIsNan<IssmDouble>(spcvx[v1]) && !xIsNan<IssmDouble>(spcvx[v2])){
-                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,
-                                                                        1,(spcvx[v1]+spcvx[v2])/(2.),DiagnosticHorizAnalysisEnum));
-                                        count++;
-                                }
-                                if(!xIsNan<IssmDouble>(spcvy[v1]) && !xIsNan<IssmDouble>(spcvy[v2])){
-                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,
-                                                                        2,(spcvy[v1]+spcvy[v2])/(2.),DiagnosticHorizAnalysisEnum));
-                                        count++;
-                                }
-                                if (reCast<int,IssmDouble>(vertices_type[v1])==FSApproximationEnum ||  (reCast<int,IssmDouble>(vertices_type[v1])==NoneApproximationEnum)){
-                                        if(!xIsNan<IssmDouble>(spcvz[v1]) && !xIsNan<IssmDouble>(spcvz[v2])){
-                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,
-                                                                                3,(spcvz[v1]+spcvz[v2])/(2.),DiagnosticHorizAnalysisEnum));
-                                                count++;
-                                        }
-                                }
-                        }
-                }
-
-                /*Clean up*/
-                xDelete<bool>(my_edges);
         }