Changeset 16343


Ignore:
Timestamp:
10/08/13 16:39:54 (11 years ago)
Author:
Mathieu Morlighem
Message:

CHG: adding some support for 2d FS

Location:
issm/trunk-jpl/src
Files:
6 added
1 deleted
22 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/Makefile.am

    r16304 r16343  
    7878                                        ./classes/Elements/ElementHook.h\
    7979                                        ./classes/Elements/ElementHook.cpp\
     80                                        ./classes/Elements/Seg.h\
     81                                        ./classes/Elements/Seg.cpp\
     82                                        ./classes/Elements/SegRef.h\
     83                                        ./classes/Elements/SegRef.cpp\
    8084                                        ./classes/Elements/Tria.h\
    8185                                        ./classes/Elements/Tria.cpp\

  • issm/trunk-jpl/src/c/analyses/stressbalance_core.cpp

    r16291 r16343  
    4747        /*Compute slopes: */
    4848        if(isSIA) surfaceslope_core(femmodel);
    49         if(isFS){
     49        if(isFS && meshtype==Mesh3DEnum){
    5050                bedslope_core(femmodel);
    5151                femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);

  • issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

    r16330 r16343  
    22422242        #endif
    22432243
    2244         //Need to know the type of approximation for this element
     2244        /*Need to know the type of approximation for this element*/
    22452245        if(iomodel->Data(FlowequationElementEquationEnum)){
    2246                 if (iomodel->Data(FlowequationElementEquationEnum)[index]==SSAApproximationEnum){
    2247                         this->inputs->AddInput(new IntInput(ApproximationEnum,SSAApproximationEnum));
    2248                 }
    2249                 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==HOApproximationEnum){
    2250                         this->inputs->AddInput(new IntInput(ApproximationEnum,HOApproximationEnum));
    2251                 }
    2252                 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==SSAHOApproximationEnum){
    2253                         this->inputs->AddInput(new IntInput(ApproximationEnum,SSAHOApproximationEnum));
    2254                 }
    2255                 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==SIAApproximationEnum){
    2256                         this->inputs->AddInput(new IntInput(ApproximationEnum,SIAApproximationEnum));
    2257                 }
    2258                 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==L1L2ApproximationEnum){
    2259                         this->inputs->AddInput(new IntInput(ApproximationEnum,L1L2ApproximationEnum));
    2260                 }
    2261                 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==FSApproximationEnum){
    2262                         this->inputs->AddInput(new IntInput(ApproximationEnum,FSApproximationEnum));
    2263                 }
    2264                 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==SSAFSApproximationEnum){
    2265                         this->inputs->AddInput(new IntInput(ApproximationEnum,SSAFSApproximationEnum));
    2266                 }
    2267                 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==HOFSApproximationEnum){
    2268                         this->inputs->AddInput(new IntInput(ApproximationEnum,HOFSApproximationEnum));
    2269                 }
    2270                 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==NoneApproximationEnum){
    2271                         this->inputs->AddInput(new IntInput(ApproximationEnum,NoneApproximationEnum));
    2272                 }
    2273                 else{
    2274                         _error_("Approximation type " << EnumToStringx(reCast<int,IssmDouble>(iomodel->Data(FlowequationElementEquationEnum)[index])) << " not supported yet");
    2275                 }
     2246                this->inputs->AddInput(new IntInput(ApproximationEnum,iomodel->Data(FlowequationElementEquationEnum)[index]));
    22762247        }
    22772248

  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r16330 r16343  
    833833}
    834834/*}}}*/
     835/*FUNCTION Tria::GetDofListVelocity{{{*/
     836void  Tria::GetDofListVelocity(int** pdoflist,int setenum){
     837
     838        /*Fetch number of nodes and dof for this finite element*/
     839        int numnodes = this->NumberofNodesVelocity();
     840
     841        /*First, figure out size of doflist and create it: */
     842        int numberofdofs=0;
     843        for(int i=0;i<numnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(FSvelocityEnum,setenum);
     844
     845        /*Allocate output*/
     846        int* doflist=xNew<int>(numberofdofs);
     847
     848        /*Populate: */
     849        int count=0;
     850        for(int i=0;i<numnodes;i++){
     851                nodes[i]->GetDofList(doflist+count,FSvelocityEnum,setenum);
     852                count+=nodes[i]->GetNumberOfDofs(FSvelocityEnum,setenum);
     853        }
     854
     855        /*Assign output pointers:*/
     856        *pdoflist=doflist;
     857}
     858/*}}}*/
     859/*FUNCTION Tria::GetDofListPressure{{{*/
     860void  Tria::GetDofListPressure(int** pdoflist,int setenum){
     861
     862        /*Fetch number of nodes and dof for this finite element*/
     863        int vnumnodes = this->NumberofNodesVelocity();
     864        int pnumnodes = this->NumberofNodesPressure();
     865
     866        /*First, figure out size of doflist and create it: */
     867        int numberofdofs=0;
     868        for(int i=vnumnodes;i<vnumnodes+pnumnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(FSApproximationEnum,setenum);
     869
     870        /*Allocate output*/
     871        int* doflist=xNew<int>(numberofdofs);
     872
     873        /*Populate: */
     874        int count=0;
     875        for(int i=vnumnodes;i<vnumnodes+pnumnodes;i++){
     876                nodes[i]->GetDofList(doflist+count,FSApproximationEnum,setenum);
     877                count+=nodes[i]->GetNumberOfDofs(FSApproximationEnum,setenum);
     878        }
     879
     880        /*Assign output pointers:*/
     881        *pdoflist=doflist;
     882}
     883/*}}}*/
    835884/*FUNCTION Tria::GetElementType {{{*/
    836885int Tria::GetElementType(){
    … …  
    13041353        #ifdef _HAVE_STRESSBALANCE_
    13051354        case StressbalanceAnalysisEnum:
    1306                 GetSolutionFromInputsStressbalanceHoriz(solution);
     1355                int approximation;
     1356                inputs->GetInputValue(&approximation,ApproximationEnum);
     1357                if(approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
     1358                        GetSolutionFromInputsStressbalanceFS(solution);
     1359                }
     1360                else if (approximation==SSAApproximationEnum){
     1361                        GetSolutionFromInputsStressbalanceHoriz(solution);
     1362                }
     1363                else{
     1364                        _error_("approximation not supported yet");
     1365                }
    13071366                break;
    13081367        case StressbalanceSIAAnalysisEnum:
    … …  
    15321591                tria_vertex_ids[i]=reCast<int>(iomodel->elements[3*index+i]); //ids for vertices are in the elements array from Matlab
    15331592        }
     1593
     1594        /*Need to know the type of approximation for this element*/
     1595        if(iomodel->Data(FlowequationElementEquationEnum)){
     1596                this->inputs->AddInput(new IntInput(ApproximationEnum,iomodel->Data(FlowequationElementEquationEnum)[index]));
     1597        }
     1598
    15341599
    15351600        /*Control Inputs*/
    … …  
    19502015                                name==VxEnum ||
    19512016                                name==VyEnum ||
     2017                                name==PressureEnum ||
    19522018                                name==InversionVxObsEnum ||
    19532019                                name==InversionVyObsEnum ||
    … …  
    24592525                        tria_node_ids[5]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+2]+1;
    24602526                        break;
     2527                case P1P1Enum: case P1P1GLSEnum:
     2528                        numnodes        = 6;
     2529                        tria_node_ids   = xNew<int>(numnodes);
     2530                        tria_node_ids[0]=iomodel->nodecounter+iomodel->elements[3*index+0];
     2531                        tria_node_ids[1]=iomodel->nodecounter+iomodel->elements[3*index+1];
     2532                        tria_node_ids[2]=iomodel->nodecounter+iomodel->elements[3*index+2];
     2533
     2534                        tria_node_ids[3]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elements[3*index+0];
     2535                        tria_node_ids[4]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elements[3*index+1];
     2536                        tria_node_ids[5]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elements[3*index+2];
     2537                        break;
     2538                case MINIEnum: case MINIcondensedEnum:
     2539                        numnodes       = 7;
     2540                        tria_node_ids  = xNew<int>(numnodes);
     2541                        tria_node_ids[0]=iomodel->nodecounter+iomodel->elements[3*index+0];
     2542                        tria_node_ids[1]=iomodel->nodecounter+iomodel->elements[3*index+1];
     2543                        tria_node_ids[2]=iomodel->nodecounter+iomodel->elements[3*index+2];
     2544                        tria_node_ids[3]=iomodel->nodecounter+iomodel->numberofvertices+index+1;
     2545
     2546                        tria_node_ids[4]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[3*index+0];
     2547                        tria_node_ids[5]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[3*index+1];
     2548                        tria_node_ids[6]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[3*index+2];
     2549                        break;
     2550                case TaylorHoodEnum:
     2551                        numnodes        = 9;
     2552                        tria_node_ids   = xNew<int>(numnodes);
     2553                        tria_node_ids[0]=iomodel->nodecounter+iomodel->elements[3*index+0];
     2554                        tria_node_ids[1]=iomodel->nodecounter+iomodel->elements[3*index+1];
     2555                        tria_node_ids[2]=iomodel->nodecounter+iomodel->elements[3*index+2];
     2556                        tria_node_ids[3]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+0]+1;
     2557                        tria_node_ids[4]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+1]+1;
     2558                        tria_node_ids[5]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+2]+1;
     2559
     2560                        tria_node_ids[6]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[3*index+0];
     2561                        tria_node_ids[7]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[3*index+1];
     2562                        tria_node_ids[8]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[3*index+2];
     2563                        break;
    24612564                default:
    24622565                        _error_("Finite element "<<EnumToStringx(finiteelement_type)<<" not supported yet");
    … …  
    33863489        delete gauss;
    33873490        return Ke;
     3491}
     3492/*}}}*/
     3493/*FUNCTION Tria::GetSolutionFromInputsStressbalanceFS{{{*/
     3494void  Tria::GetSolutionFromInputsStressbalanceFS(Vector<IssmDouble>* solution){
     3495
     3496        int*         vdoflist=NULL;
     3497        int*         pdoflist=NULL;
     3498        IssmDouble   vx,vy,p;
     3499        IssmDouble   FSreconditioning;
     3500        GaussTria   *gauss;
     3501
     3502        /*Fetch number of nodes and dof for this finite element*/
     3503        int vnumnodes = this->NumberofNodesVelocity();
     3504        int pnumnodes = this->NumberofNodesPressure();
     3505        int vnumdof   = vnumnodes*NDOF2;
     3506        int pnumdof   = pnumnodes*NDOF1;
     3507
     3508        /*Initialize values*/
     3509        IssmDouble* vvalues = xNew<IssmDouble>(vnumdof);
     3510        IssmDouble* pvalues = xNew<IssmDouble>(pnumdof);
     3511
     3512        /*Get dof list: */
     3513        GetDofListVelocity(&vdoflist,GsetEnum);
     3514        GetDofListPressure(&pdoflist,GsetEnum);
     3515        Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
     3516        Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
     3517        Input* p_input =inputs->GetInput(PressureEnum); _assert_(p_input);
     3518
     3519        this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
     3520
     3521        /*Ok, we have vx vy vz in values, fill in vx vy vz arrays: */
     3522        gauss = new GaussTria();
     3523        for(int i=0;i<vnumnodes;i++){
     3524                gauss->GaussNode(this->VelocityInterpolation(),i);
     3525                vx_input->GetInputValue(&vx,gauss);
     3526                vy_input->GetInputValue(&vy,gauss);
     3527                vvalues[i*NDOF2+0]=vx;
     3528                vvalues[i*NDOF2+1]=vy;
     3529        }
     3530        for(int i=0;i<pnumnodes;i++){
     3531                gauss->GaussNode(this->PressureInterpolation(),i);
     3532                p_input->GetInputValue(&p ,gauss);
     3533                pvalues[i]=p/FSreconditioning;
     3534        }
     3535
     3536        /*Add value to global vector*/
     3537        solution->SetValues(vnumdof,vdoflist,vvalues,INS_VAL);
     3538        solution->SetValues(pnumdof,pdoflist,pvalues,INS_VAL);
     3539
     3540        /*Free ressources:*/
     3541        delete gauss;
     3542        xDelete<int>(pdoflist);
     3543        xDelete<int>(vdoflist);
     3544        xDelete<IssmDouble>(pvalues);
     3545        xDelete<IssmDouble>(vvalues);
    33883546}
    33893547/*}}}*/

  • issm/trunk-jpl/src/c/classes/Elements/Tria.h

    r16314 r16343  
    205205                int            GetElementType(void);
    206206                void             GetDofList(int** pdoflist,int approximation_enum,int setenum);
     207                void             GetDofListVelocity(int** pdoflist,int setenum);
     208                void             GetDofListPressure(int** pdoflist,int setenum);
    207209                void             GetVertexPidList(int* doflist);
    208210                void           GetVertexSidList(int* sidlist);
    … …  
    236238                ElementVector* CreatePVectorStressbalanceSIA(void);
    237239                ElementMatrix* CreateJacobianStressbalanceSSA(void);
     240                void      GetSolutionFromInputsStressbalanceFS(Vector<IssmDouble>* solution);
    238241                void      GetSolutionFromInputsStressbalanceHoriz(Vector<IssmDouble>* solution);
    239242                void      GetSolutionFromInputsStressbalanceSIA(Vector<IssmDouble>* solution);

  • issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp

    r15767 r16343  
    4343/*FUNCTION TriaRef::SetElementType{{{*/
    4444void TriaRef::SetElementType(int type,int type_counter){
    45 
    46         _assert_(type==P1Enum || type==P1DGEnum || type==P1bubbleEnum || type==P1bubblecondensedEnum || type==P2Enum);
    4745
    4846        /*initialize element type*/
    … …  
    653651                case P1bubblecondensedEnum: return NUMNODESP1b;
    654652                case P2Enum:                return NUMNODESP2;
     653                case P1P1Enum:              return NUMNODESP1*2;
     654                case P1P1GLSEnum:           return NUMNODESP1*2;
     655                case MINIcondensedEnum:     return NUMNODESP1b+NUMNODESP1;
     656                case MINIEnum:              return NUMNODESP1b+NUMNODESP1;
     657                case TaylorHoodEnum:        return NUMNODESP2+NUMNODESP1;
    655658                default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
    656659        }
    … …  
    659662}
    660663/*}}}*/
     664/*FUNCTION TriaRef::NumberofNodesPressure{{{*/
     665int TriaRef::NumberofNodesPressure(void){
     666
     667        switch(this->element_type){
     668                case P1P1Enum:          return NUMNODESP1;
     669                case P1P1GLSEnum:       return NUMNODESP1;
     670                case MINIcondensedEnum: return NUMNODESP1;
     671                case MINIEnum:          return NUMNODESP1;
     672                case TaylorHoodEnum:    return NUMNODESP1;
     673                default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     674        }
     675
     676        return -1;
     677}
     678/*}}}*/
     679/*FUNCTION TriaRef::NumberofNodesVelocity{{{*/
     680int TriaRef::NumberofNodesVelocity(void){
     681
     682        switch(this->element_type){
     683                case P1P1Enum:          return NUMNODESP1;
     684                case P1P1GLSEnum:       return NUMNODESP1;
     685                case MINIcondensedEnum: return NUMNODESP1b;
     686                case MINIEnum:          return NUMNODESP1b;
     687                case TaylorHoodEnum:    return NUMNODESP2;
     688                default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     689        }
     690
     691        return -1;
     692}
     693/*}}}*/
     694/*FUNCTION TriaRef::VelocityInterpolation{{{*/
     695int TriaRef::VelocityInterpolation(void){
     696
     697        switch(this->element_type){
     698                case P1P1Enum:          return P1Enum;
     699                case P1P1GLSEnum:       return P1Enum;
     700                case MINIcondensedEnum: return P1bubbleEnum;
     701                case MINIEnum:          return P1bubbleEnum;
     702                case TaylorHoodEnum:    return P2Enum;
     703                default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     704        }
     705
     706        return -1;
     707}
     708/*}}}*/
     709/*FUNCTION TriaRef::PressureInterpolation{{{*/
     710int TriaRef::PressureInterpolation(void){
     711
     712        switch(this->element_type){
     713                case P1P1Enum:          return P1Enum;
     714                case P1P1GLSEnum:       return P1Enum;
     715                case MINIcondensedEnum: return P1Enum;
     716                case MINIEnum:          return P1Enum;
     717                case TaylorHoodEnum:    return P1Enum;
     718                default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     719        }
     720
     721        return -1;
     722}
     723/*}}}*/

  • issm/trunk-jpl/src/c/classes/Elements/TriaRef.h

    r15767 r16343  
    4545
    4646                int  NumberofNodes(void);
     47                int  NumberofNodesVelocity(void);
     48                int  NumberofNodesPressure(void);
     49                int  VelocityInterpolation(void);
     50                int  PressureInterpolation(void);
    4751};
    4852#endif

  • issm/trunk-jpl/src/c/classes/Node.cpp

    r16338 r16343  
    3939        /*indexing:*/
    4040        this->indexingupdate = true;
    41         DistributeNumDofs(&this->indexing,analysis_type,in_approximation); //number of dofs per node
     41        DistributeNumDofs(&this->indexing,analysis_type,in_approximation,iomodel->meshtype); //number of dofs per node
    4242
    4343        if(analysis_type==StressbalanceAnalysisEnum)
    … …  
    920920                        case PressureEnum: numdofs+=1; break;
    921921                        case XYEnum:       numdofs+=2; break;
    922                         case XZEnum:       numdofs+=2; break;
    923922                        case XYZEnum:      numdofs+=3; break;
    924923                        default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
    … …  
    971970                        case PressureEnum: numdofs+=1; break;
    972971                        case XYEnum:       numdofs+=2; break;
    973                         case XZEnum:       numdofs+=2; break;
    974972                        case XYZEnum:      numdofs+=3; break;
    975973                        default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
    … …  
    10211019                        case PressureEnum: numdofs+=1; break;
    10221020                        case XYEnum:       numdofs+=2; break;
    1023                         case XZEnum:       numdofs+=2; break;
    10241021                        case XYZEnum:      numdofs+=3; break;
    10251022                        default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
    … …  
    10701067                        case PressureEnum: numdofs+=1; break;
    10711068                        case XYEnum:       numdofs+=2; break;
    1072                         case XZEnum:       numdofs+=2; break;
    10731069                        case XYZEnum:      numdofs+=3; break;
    10741070                        default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
    … …  
    11101106                        case PressureEnum: numdofs+=1; break;
    11111107                        case XYEnum:       numdofs+=2; break;
    1112                         case XZEnum:       numdofs+=2; break;
    11131108                        case XYZEnum:      numdofs+=3; break;
    11141109                        default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
    … …  
    11471142                                counter+=2;
    11481143                                break;
    1149                         case XZEnum:
    1150                                 /*We remove the y component, we need to renormalize x and z: x=[x1 0 x2] y=[-x2 0 x1]*/
    1151                                 norm = sqrt( coord_system[0][0]*coord_system[0][0] + coord_system[2][0]*coord_system[2][0]); _assert_(norm>1.e-4);
    1152                                 transform[(numdofs)*(counter+0) + counter+0] =   coord_system[0][0]/norm;
    1153                                 transform[(numdofs)*(counter+0) + counter+2] = - coord_system[2][0]/norm;
    1154                                 transform[(numdofs)*(counter+1) + counter+0] =   coord_system[2][0]/norm;
    1155                                 transform[(numdofs)*(counter+1) + counter+2] =   coord_system[0][0]/norm;
    1156                                 counter+=2;
    1157                                 break;
    11581144                        case XYZEnum:
    11591145                                /*The 3 coordinates are changed (x,y,z)*/

  • issm/trunk-jpl/src/c/classes/classes.h

    r16304 r16343  
    3939#include "./Elements/Penta.h"
    4040#include "./Elements/PentaRef.h"
     41#include "./Elements/Seg.h"
     42#include "./Elements/SegRef.h"
    4143#include "./Elements/Tria.h"
    4244#include "./Elements/TriaRef.h"

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp

    r16181 r16343  
    66#include "../../classes/classes.h"
    77
    8 void DistributeNumDofs(DofIndexing* index,int analysis_type,int node_type){
     8void DistributeNumDofs(DofIndexing* index,int analysis_type,int node_type,int mesh_type){
    99
    1010        /*For now, we distribute by analysis_type, later, we will distribute using the analysis_type,
    … …  
    3131                                        break;
    3232                                case FSvelocityEnum:
    33                                         numdofs=3;
     33                                        if(mesh_type==Mesh3DEnum){
     34                                                numdofs=3;
     35                                        }
     36                                        else if(mesh_type==Mesh2DverticalEnum){
     37                                                numdofs=2;
     38                                        }
     39                                        else{
     40                                                _error_("mesh type not supported yet");
     41                                        }
    3442                                        break;
    3543                                case FSpressureEnum:
    … …  
    3745                                        break;
    3846                                case NoneApproximationEnum:
    39                                         numdofs=4;
     47                                        if(mesh_type==Mesh3DEnum){
     48                                                numdofs=4;
     49                                        }
     50                                        else if(mesh_type==Mesh2DverticalEnum){
     51                                                numdofs=3;
     52                                        }
     53                                        else{
     54                                                _error_("mesh type not supported yet");
     55                                        }
    4056                                        break;
    4157                                case SSAHOApproximationEnum:

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h

    r16304 r16343  
    158158
    159159/*Distribution of dofs: */
    160 void DistributeNumDofs(DofIndexing* index,int analysis_type,int node_type);
     160void DistributeNumDofs(DofIndexing* index,int analysis_type,int node_type,int mesh_type);
    161161
    162162#endif

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/Stressbalance/CreateConstraintsStressbalance.cpp

    r16291 r16343  
    8888                }
    8989
    90                 IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,1);
    91                 IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvyEnum,StressbalanceAnalysisEnum,finiteelement,2);
    92 
    9390                if(isFS){
    9491
    9592                        /*Constraint at the bedrock interface (v.n = vz = 0) (Coordinates will be updated according to the bed slope)*/
    96                         iomodel->FetchData(&spcvz,&Mz,&Nz,StressbalanceSpcvzEnum);
    9793                        iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum);
    9894                        iomodel->FetchData(&nodeonFS,NULL,NULL,FlowequationBorderFSEnum);
    9995                        iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum);
    10096                        iomodel->FetchData(&groundedice_ls,NULL,NULL,MaskGroundediceLevelsetEnum);
     97                        if(iomodel->meshtype==Mesh3DEnum){
     98                                iomodel->FetchData(&spcvz,&Mz,&Nz,StressbalanceSpcvzEnum);
     99                        }
     100                        else if (iomodel->meshtype==Mesh2DverticalEnum){
     101                                iomodel->FetchData(&spcvz,&Mz,&Nz,StressbalanceSpcvyEnum);
     102                        }
     103                        else{
     104                                _error_("not supported yet");
     105                        }
    101106                        for(i=0;i<iomodel->numberofvertices;i++){
    102107                                if(iomodel->my_vertices[i]){
    … …  
    111116                                }
    112117                        }
    113                         IoModelToConstraintsx(constraints,iomodel,spcvz,Mz,Nz,StressbalanceAnalysisEnum,finiteelement,3);
    114                         iomodel->DeleteData(spcvz,StressbalanceSpcvzEnum);
     118                        if(iomodel->meshtype==Mesh3DEnum){
     119                                IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,1);
     120                                IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvyEnum,StressbalanceAnalysisEnum,finiteelement,2);
     121                                IoModelToConstraintsx(constraints,iomodel,spcvz,Mz,Nz,StressbalanceAnalysisEnum,finiteelement,3);
     122                                iomodel->DeleteData(spcvz,StressbalanceSpcvzEnum);
     123                        }
     124                        else if (iomodel->meshtype==Mesh2DverticalEnum){
     125                                IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,1);
     126                                IoModelToConstraintsx(constraints,iomodel,spcvz,Mz,Nz,StressbalanceAnalysisEnum,finiteelement,2);
     127                                iomodel->DeleteData(spcvz,StressbalanceSpcvyEnum);
     128                        }
     129                        else{
     130                                _error_("not supported yet");
     131                        }
    115132                        iomodel->DeleteData(vertices_type,FlowequationVertexEquationEnum);
    116133                        iomodel->DeleteData(nodeonFS,FlowequationBorderFSEnum);
    … …  
    151168                        iomodel->DeleteData(z,MeshZEnum);
    152169                }
     170                else{
     171                        IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,1);
     172                        IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvyEnum,StressbalanceAnalysisEnum,finiteelement,2);
     173                }
    153174
    154175                *pconstraints=constraints;

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/Stressbalance/UpdateElementsStressbalance.cpp

    r16291 r16343  
    112112                iomodel->FetchDataToInput(elements,VzEnum,0.);
    113113                if(dakota_analysis)elements->InputDuplicate(VzEnum,QmuVzEnum);
    114                 if(isFS){
    115                         iomodel->FetchDataToInput(elements,PressureEnum,0.);
    116                         if(dakota_analysis)elements->InputDuplicate(PressureEnum,QmuPressureEnum);
    117                 }
     114        }
     115        if(isFS){
     116                iomodel->FetchDataToInput(elements,PressureEnum,0.);
     117                if(dakota_analysis)elements->InputDuplicate(PressureEnum,QmuPressureEnum);
    118118        }
    119119

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/StressbalanceVertical/CreateConstraintsStressbalanceVertical.cpp

    r16291 r16343  
    2020        Constraints* constraints=*pconstraints;
    2121
    22         /*return if 2d mesh*/
    23         if(iomodel->meshtype==Mesh2DhorizontalEnum) return;
     22        /*return if not 3d mesh*/
     23        if(iomodel->meshtype!=Mesh3DEnum) return;
    2424
    2525        /*Fetch data: */

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/StressbalanceVertical/CreateNodesStressbalanceVertical.cpp

    r16291 r16343  
    1111void    CreateNodesStressbalanceVertical(Nodes** pnodes, IoModel* iomodel){
    1212
    13         /*Now, is the flag macayaealHO on? otherwise, do nothing: */
    14         if(iomodel->meshtype==Mesh2DhorizontalEnum) return;
     13        /*return if not 3d mesh*/
     14        if(iomodel->meshtype!=Mesh3DEnum) return;
    1515
    1616        iomodel->FetchData(3,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationVertexEquationEnum);

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/StressbalanceVertical/UpdateElementsStressbalanceVertical.cpp

    r16291 r16343  
    1111void    UpdateElementsStressbalanceVertical(Elements* elements, IoModel* iomodel,int analysis_counter,int analysis_type){
    1212
    13         /*Now, is the model 3d? otherwise, do nothing: */
    14         if (iomodel->meshtype==Mesh2DhorizontalEnum)return;
     13        /*return if not 3d mesh*/
     14        if(iomodel->meshtype!=Mesh3DEnum) return;
    1515
    1616        /*Update elements: */

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r16338 r16343  
    403403        StringArrayParamEnum,
    404404        StringParamEnum,
     405        SegEnum,
     406        SegInputEnum,
    405407        TriaEnum,
    406408        TriaInputEnum,
    … …  
    611613        /*Coordinate Systems{{{*/
    612614        XYEnum,
    613         XZEnum,
    614615        XYZEnum,
    615616        /*}}}*/

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r16338 r16343  
    403403                case StringArrayParamEnum : return "StringArrayParam";
    404404                case StringParamEnum : return "StringParam";
     405                case SegEnum : return "Seg";
     406                case SegInputEnum : return "SegInput";
    405407                case TriaEnum : return "Tria";
    406408                case TriaInputEnum : return "TriaInput";
    … …  
    587589                case NearestInterpEnum : return "NearestInterp";
    588590                case XYEnum : return "XY";
    589                 case XZEnum : return "XZ";
    590591                case XYZEnum : return "XYZ";
    591592                case DenseEnum : return "Dense";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r16338 r16343  
    412412              else if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum;
    413413              else if (strcmp(name,"StringParam")==0) return StringParamEnum;
     414              else if (strcmp(name,"Seg")==0) return SegEnum;
     415              else if (strcmp(name,"SegInput")==0) return SegInputEnum;
    414416              else if (strcmp(name,"Tria")==0) return TriaEnum;
    415417              else if (strcmp(name,"TriaInput")==0) return TriaInputEnum;
    … …  
    504506              else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum;
    505507              else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum;
    506               else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
    507               else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum;
     511              if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
     512              else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum;
     513              else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum;
    512514              else if (strcmp(name,"GiaCrossSectionShape")==0) return GiaCrossSectionShapeEnum;
    513515              else if (strcmp(name,"GiadWdt")==0) return GiadWdtEnum;
    … …  
    599601              else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum;
    600602              else if (strcmp(name,"XY")==0) return XYEnum;
    601               else if (strcmp(name,"XZ")==0) return XZEnum;
    602603              else if (strcmp(name,"XYZ")==0) return XYZEnum;
    603604              else if (strcmp(name,"Dense")==0) return DenseEnum;

  • issm/trunk-jpl/src/m/classes/mesh2dvertical.m

    r16335 r16343  
    77        properties (SetAccess=public)
    88                x                           = NaN;
    9                 z                           = NaN;
     9                y                           = NaN;
    1010                elements                    = NaN
    1111                numberofelements            = 0;
    … …  
    4848
    4949                        md = checkfield(md,'mesh.x','NaN',1,'size',[md.mesh.numberofvertices 1]);
    50                         md = checkfield(md,'mesh.z','NaN',1,'size',[md.mesh.numberofvertices 1]);
     50                        md = checkfield(md,'mesh.y','NaN',1,'size',[md.mesh.numberofvertices 1]);
    5151                        md = checkfield(md,'mesh.elements','NaN',1,'>',0,'values',1:md.mesh.numberofvertices);
    5252                        md = checkfield(md,'mesh.elements','size',[md.mesh.numberofelements 3]);
    … …  
    7272                        fielddisplay(obj,'elements','vertex indices of the mesh elements');
    7373                        fielddisplay(obj,'x','vertices x coordinate [m]');
    74                         fielddisplay(obj,'z','vertices z coordinate [m]');
     74                        fielddisplay(obj,'y','vertices y coordinate [m]');
    7575                        fielddisplay(obj,'edges','edges of the 2d mesh (vertex1 vertex2 element1 element2)');
    7676                        fielddisplay(obj,'numberofedges','number of edges of the 2d mesh');
    … …  
    8585                        fielddisplay(obj,'average_vertex_connectivity','average number of vertices connected to one vertex');
    8686
    87                         disp(sprintf('\n      Extracted model:'));
    88                         fielddisplay(obj,'extractedvertices','vertices extracted from the model');
    89                         fielddisplay(obj,'extractedelements','elements extracted from the model');
    90 
    9187                        disp(sprintf('\n      Projection:'));
    9288                        fielddisplay(obj,'lat','vertices latitude [degrees]');
    … …  
    9793                        WriteData(fid,'enum',MeshTypeEnum(),'data',StringToEnum(['Mesh' meshtype(obj)]),'format','Integer');
    9894                        WriteData(fid,'object',obj,'class','mesh','fieldname','x','format','DoubleMat','mattype',1);
    99                         WriteData(fid,'object',obj,'class','mesh','fieldname','z','format','DoubleMat','mattype',1);
    100                         WriteData(fid,'enum',MeshYEnum(),'data',zeros(obj.numberofvertices,1),'format','DoubleMat','mattype',1);
     95                        WriteData(fid,'object',obj,'class','mesh','fieldname','y','format','DoubleMat','mattype',1);
     96                        WriteData(fid,'enum',MeshZEnum(),'data',zeros(obj.numberofvertices,1),'format','DoubleMat','mattype',1);
    10197                        WriteData(fid,'object',obj,'class','mesh','fieldname','elements','format','DoubleMat','mattype',2);
    10298                        WriteData(fid,'object',obj,'class','mesh','fieldname','numberofelements','format','Integer');
    … …  
    207203                        elements = self.elements;
    208204                        x        = self.x;
    209                         y        = self.z;
     205                        y        = self.y;
    210206                        z        = zeros(self.numberofvertices,1);
    211207                end % }}}

  • issm/trunk-jpl/src/m/enum/EnumDefinitions.py

    r16338 r16343  
    395395def StringArrayParamEnum(): return StringToEnum("StringArrayParam")[0]
    396396def StringParamEnum(): return StringToEnum("StringParam")[0]
     397def SegEnum(): return StringToEnum("Seg")[0]
     398def SegInputEnum(): return StringToEnum("SegInput")[0]
    397399def TriaEnum(): return StringToEnum("Tria")[0]
    398400def TriaInputEnum(): return StringToEnum("TriaInput")[0]
    … …  
    579581def NearestInterpEnum(): return StringToEnum("NearestInterp")[0]
    580582def XYEnum(): return StringToEnum("XY")[0]
    581 def XZEnum(): return StringToEnum("XZ")[0]
    582583def XYZEnum(): return StringToEnum("XYZ")[0]
    583584def DenseEnum(): return StringToEnum("Dense")[0]

  • issm/trunk-jpl/src/m/mesh/bamg.m

    r16335 r16343  
    328328        md.mesh=mesh2dvertical;
    329329        md.mesh.x=bamgmesh_out.Vertices(:,1);
    330         md.mesh.z=bamgmesh_out.Vertices(:,2);
     330        md.mesh.y=bamgmesh_out.Vertices(:,2);
    331331        md.mesh.elements=bamgmesh_out.Triangles(:,1:3);
    332332        md.mesh.edges=bamgmesh_out.IssmEdges;
Note: See TracChangeset for help on using the changeset viewer.