Changeset 15636


Ignore:
Timestamp:
07/26/13 09:24:48 (12 years ago)
Author:
Mathieu Morlighem
Message:

NEW: added P1xP2 element for Penta, needs to be debugged

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

Legend:

Unmodified
Added
Removed

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

    r15634 r15636  
    31773177                        penta_node_ids[4]=iomodel->nodecounter+iomodel->elements[6*index+4];
    31783178                        penta_node_ids[5]=iomodel->nodecounter+iomodel->elements[6*index+5];
     3179                        break;
     3180                case P1xP2Enum:
     3181                        numnodes         = 9;
     3182                        penta_node_ids   = xNew<int>(numnodes);
     3183                        penta_node_ids[ 0]=iomodel->nodecounter+iomodel->elements[6*index+0];
     3184                        penta_node_ids[ 1]=iomodel->nodecounter+iomodel->elements[6*index+1];
     3185                        penta_node_ids[ 2]=iomodel->nodecounter+iomodel->elements[6*index+2];
     3186                        penta_node_ids[ 3]=iomodel->nodecounter+iomodel->elements[6*index+3];
     3187                        penta_node_ids[ 4]=iomodel->nodecounter+iomodel->elements[6*index+4];
     3188                        penta_node_ids[ 5]=iomodel->nodecounter+iomodel->elements[6*index+5];
     3189                        penta_node_ids[ 6]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+0]+1;
     3190                        penta_node_ids[ 7]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+1]+1;
     3191                        penta_node_ids[ 8]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+2]+1;
    31793192                        break;
    31803193                case P2xP1Enum:

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

    r15627 r15636  
    1919#define NUMNODESP1_2d 3
    2020#define NUMNODESP1b   7
     21#define NUMNODESP1xP2 9
    2122#define NUMNODESP2xP1 12
    2223#define NUMNODESP2    15
    … …  
    11291130                        basis[10]=4.*gauss->coord3*gauss->coord1*(1.+gauss->coord4)/2.;
    11301131                        basis[11]=4.*gauss->coord1*gauss->coord2*(1.+gauss->coord4)/2.;
     1132                        return;
     1133                case P1xP2Enum:
     1134                        /*Corner nodes*/
     1135                        basis[ 0]=gauss->coord1*gauss->coord4*(gauss->coord4-1.)/2.;
     1136                        basis[ 1]=gauss->coord2*gauss->coord4*(gauss->coord4-1.)/2.;
     1137                        basis[ 2]=gauss->coord3*gauss->coord4*(gauss->coord4-1.)/2.;
     1138                        basis[ 3]=gauss->coord1*gauss->coord4*(gauss->coord4+1.)/2.;
     1139                        basis[ 4]=gauss->coord2*gauss->coord4*(gauss->coord4+1.)/2.;
     1140                        basis[ 5]=gauss->coord3*gauss->coord4*(gauss->coord4+1.)/2.;
     1141                        /*mid-sides of quads*/
     1142                        basis[ 6]=gauss->coord1*(1.-gauss->coord4*gauss->coord4);
     1143                        basis[ 7]=gauss->coord2*(1.-gauss->coord4*gauss->coord4);
     1144                        basis[ 8]=gauss->coord3*(1.-gauss->coord4*gauss->coord4);
    11311145                        return;
    11321146                case P2Enum:
    … …  
    13091323                        dbasis[NUMNODESP2xP1*2+11] = 2.*gauss->coord1*gauss->coord2;
    13101324                        return;
     1325                case P1xP2Enum:
     1326                        /*Nodal function 1*/
     1327                        dbasis[NUMNODESP1xP2*0+0]   = zeta*(zeta-1.)/4.;
     1328                        dbasis[NUMNODESP1xP2*1+0]   = SQRT3/12.*zeta*(zeta-1.);
     1329                        dbasis[NUMNODESP1xP2*2+0]   = .5*(1.-2.*zeta)*gauss->coord1;
     1330                        /*Nodal function 2*/
     1331                        dbasis[NUMNODESP1xP2*0+1]   = zeta*(1.-zeta)/4.;
     1332                        dbasis[NUMNODESP1xP2*1+1]   = SQRT3/12.*zeta*(zeta-1);
     1333                        dbasis[NUMNODESP1xP2*2+1]   = .5*(1.-2.*zeta)*gauss->coord2;
     1334                        /*Nodal function 3*/
     1335                        dbasis[NUMNODESP1xP2*0+2]   = 0.;
     1336                        dbasis[NUMNODESP1xP2*1+2]   = SQRT3/6.*zeta*(1.-zeta);
     1337                        dbasis[NUMNODESP1xP2*2+2]   = .5*(1.-2.*zeta)*gauss->coord3;
     1338                        /*Nodal function 4*/
     1339                        dbasis[NUMNODESP1xP2*0+3]   = -zeta*(1.+zeta)/4.;
     1340                        dbasis[NUMNODESP1xP2*1+3]   = -SQRT3/12.*(1.+zeta);
     1341                        dbasis[NUMNODESP1xP2*2+3]   = .5*(1.+2.*zeta)*gauss->coord1;
     1342                        /*Nodal function 5*/
     1343                        dbasis[NUMNODESP1xP2*0+4]   = zeta*(1.+zeta)/4.;
     1344                        dbasis[NUMNODESP1xP2*1+4]   = -SQRT3/12.*zeta*(1.+zeta);
     1345                        dbasis[NUMNODESP1xP2*2+4]   = .5*(1.+2.*zeta)*gauss->coord2;
     1346                        /*Nodal function 6*/
     1347                        dbasis[NUMNODESP1xP2*0+5]   = 0.;
     1348                        dbasis[NUMNODESP1xP2*1+5]   = SQRT3/6.*zeta*(1.+zeta);
     1349                        dbasis[NUMNODESP1xP2*2+5]   = .5*(1.+2.*zeta)*gauss->coord3;
     1350
     1351                        /*Nodal function 7*/
     1352                        dbasis[NUMNODESP1xP2*0+6 ] = -0.5*(1.-zeta*zeta);
     1353                        dbasis[NUMNODESP1xP2*1+6 ] = -SQRT3/6.*(1.-zeta*zeta);
     1354                        dbasis[NUMNODESP1xP2*2+6 ] = -2.*zeta*gauss->coord1;
     1355                        /*Nodal function 8*/
     1356                        dbasis[NUMNODESP1xP2*0+7 ] = 0.5*(1.-zeta*zeta);
     1357                        dbasis[NUMNODESP1xP2*1+7 ] = -SQRT3/6.*(1.-zeta*zeta);
     1358                        dbasis[NUMNODESP1xP2*2+7 ] = -2.*zeta*gauss->coord2;
     1359                        /*Nodal function 9*/
     1360                        dbasis[NUMNODESP1xP2*0+8 ] = 0.;
     1361                        dbasis[NUMNODESP1xP2*1+8 ] = SQRT3/3.*(1.-zeta*zeta);
     1362                        dbasis[NUMNODESP1xP2*2+8 ] = -2.*zeta*gauss->coord3;
     1363                        return;
    13111364                case P2Enum:
    13121365                        /*Nodal function 1*/
    … …  
    16281681                case P2Enum:    return NUMNODESP2;
    16291682                case P2xP1Enum: return NUMNODESP2xP1;
     1683                case P1xP2Enum: return NUMNODESP1xP2;
    16301684                case MINIEnum:  return NUMNODESP1b;
    16311685                default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");

  • issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp

    r15627 r15636  
    415415                        }
    416416                        break;
     417                case 9: //P1xP2 Lagrange element
     418                        switch(iv){
     419                                case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
     420                                case 1: coord1=0.; coord2=1.; coord3=0.; coord4=-1.; break;
     421                                case 2: coord1=0.; coord2=0.; coord3=1.; coord4=-1.; break;
     422                                case 3: coord1=1.; coord2=0.; coord3=0.; coord4=+1.; break;
     423                                case 4: coord1=0.; coord2=1.; coord3=0.; coord4=+1.; break;
     424                                case 5: coord1=0.; coord2=0.; coord3=1.; coord4=+1.; break;
     425
     426                                case 6: coord1=1.; coord2=0.; coord3=0.; coord4=0.; break;
     427                                case 7: coord1=0.; coord2=1.; coord3=0.; coord4=0.; break;
     428                                case 8: coord1=0.; coord2=0.; coord3=1.; coord4=0.; break;
     429                                default: _error_("node index should be in [0 5]");
     430                        }
     431                        break;
    417432                case 12: //P2xP1 Lagrange element
    418433                        switch(iv){

  • issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp

    r15635 r15636  
    4242                        /*Nothing else to do*/
    4343                        break;
    44                 case P2Enum:
     44                case P1xP2Enum:
    4545                        EdgesPartitioning(&my_edges,iomodel);
    4646                        break;
    4747                case P2xP1Enum:
     48                        EdgesPartitioning(&my_edges,iomodel);
     49                        break;
     50                case P2Enum:
    4851                        EdgesPartitioning(&my_edges,iomodel);
    4952                        break;
    … …  
    8285                                                                                        dof,(spcdata[v1]+spcdata[v2])/2.,analysis_type));
    8386                                                        count++;
     87                                                }
     88                                        }
     89                                }
     90                                break;
     91                        case P1xP2Enum:
     92                                for(i=0;i<iomodel->numberofvertices;i++){
     93                                        if((iomodel->my_vertices[i])){
     94                                                if (!xIsNan<IssmDouble>(spcdata[i])){
     95                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,dof,spcdata[i],analysis_type));
     96                                                        count++;
     97                                                }
     98                                        }
     99                                }
     100                                for(i=0;i<iomodel->numberofedges;i++){
     101                                        if(iomodel->edges[i*3+2]==2){
     102                                                if(my_edges[i]){
     103                                                        v1 = iomodel->edges[3*i+0]-1;
     104                                                        v2 = iomodel->edges[3*i+1]-1;
     105                                                        if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
     106                                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,
     107                                                                                                dof,(spcdata[v1]+spcdata[v2])/2.,analysis_type));
     108                                                                count++;
     109                                                        }
    84110                                                }
    85111                                        }
    … …  
    177203                                                }
    178204                                                xDelete<IssmDouble>(values);
     205                                        }
     206                                }
     207                                break;
     208                        case P1xP2Enum:
     209                                for(i=0;i<iomodel->numberofvertices;i++){
     210                                        if((iomodel->my_vertices[i])){
     211
     212                                                /*figure out times and values: */
     213                                                values=xNew<IssmDouble>(N);
     214                                                spcpresent=false;
     215                                                for(j=0;j<N;j++){
     216                                                        values[j]=spcdata[i*N+j];
     217                                                        if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
     218                                                }
     219
     220                                                if(spcpresent){
     221                                                        constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,dof,N,times,values,analysis_type));
     222                                                        count++;
     223                                                }
     224                                                xDelete<IssmDouble>(values);
     225                                        }
     226                                }
     227                                for(i=0;i<iomodel->numberofedges;i++){
     228                                        if(iomodel->edges[i*3+2]==2){
     229                                                if(my_edges[i]){
     230                                                        v1 = iomodel->edges[3*i+0]-1;
     231                                                        v2 = iomodel->edges[3*i+1]-1;
     232                                                        values=xNew<IssmDouble>(N);
     233                                                        spcpresent=false;
     234                                                        for(j=0;j<N;j++){
     235                                                                values[j]=(spcdata[v1*N+j]+spcdata[v2*N+j])/2.;
     236                                                                if(!xIsNan<IssmDouble>(values[j])) spcpresent=true; //NaN means no spc by default
     237                                                        }
     238                                                        if(spcpresent){
     239                                                                constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,dof,
     240                                                                                                N,times,values,analysis_type));
     241                                                                count++;
     242                                                        }
     243                                                        xDelete<IssmDouble>(values);
     244                                                }
    179245                                        }
    180246                                }

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

    r15634 r15636  
    4444                        break;
    4545
     46                case P1xP2Enum:
     47                        EdgesPartitioning(&my_edges,iomodel);
     48                        for(i=0;i<iomodel->numberofvertices;i++){
     49                                if(iomodel->my_vertices[i]){
     50                                        nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i,iomodel,analysis,approximation));
     51                                }
     52                        }
     53
     54                        counter = iomodel->numberofvertices;
     55                        for(i=0;i<iomodel->numberofedges;i++){
     56                                if(iomodel->edges[i*3+2]==2){
     57                                        if(my_edges[i]){
     58                                                node = new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,counter+1,0,iomodel,analysis,approximation);
     59                                                nodes->AddObject(node);
     60                                        }
     61                                        counter++;
     62                                }
     63                        }
     64                        break;
    4665                case P2xP1Enum:
    4766                        EdgesPartitioning(&my_edges,iomodel);

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

    r15627 r15636  
    8888                                case 1 : finiteelement = P2Enum;    break;
    8989                                case 2 : finiteelement = P2xP1Enum; break;
     90                                case 3 : finiteelement = P1xP2Enum; break;
    9091                                default: _error_("finite element "<<temp<<" not supported");
    9192                        }

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

    r15627 r15636  
    5454                                case 1 : finiteelement = P2Enum;    break;
    5555                                case 2 : finiteelement = P2xP1Enum; break;
     56                                case 3 : finiteelement = P1xP2Enum; break;
    5657                                default: _error_("finite element "<<temp<<" not supported");
    5758                        }

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

    r15627 r15636  
    5959                                case 1 : finiteelement = P2Enum;    break;
    6060                                case 2 : finiteelement = P2xP1Enum; break;
     61                                case 3 : finiteelement = P1xP2Enum; break;
    6162                                default: _error_("finite element "<<temp<<" not supported");
    6263                        }

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

    r15621 r15636  
    503503        P2Enum,
    504504        P2xP1Enum,
     505        P1xP2Enum,
    505506        MINIEnum,
    506507        MINIcondensedEnum,

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

    r15621 r15636  
    495495                case P2Enum : return "P2";
    496496                case P2xP1Enum : return "P2xP1";
     497                case P1xP2Enum : return "P1xP2";
    497498                case MINIEnum : return "MINI";
    498499                case MINIcondensedEnum : return "MINIcondensed";

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

    r15621 r15636  
    504504              else if (strcmp(name,"P2")==0) return P2Enum;
    505505              else if (strcmp(name,"P2xP1")==0) return P2xP1Enum;
     506              else if (strcmp(name,"P1xP2")==0) return P1xP2Enum;
    506507              else if (strcmp(name,"MINI")==0) return MINIEnum;
    507               else if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum;
     511              if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum;
     512              else if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum;
    512513              else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum;
    513514              else if (strcmp(name,"BoolElementResult")==0) return BoolElementResultEnum;

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

    r15627 r15636  
    8080                                md = checkfield(md,'flowequation.isFS','numel',[1],'values',[0 1]);
    8181                                md = checkfield(md,'flowequation.fe_SSA','numel',[1],'values',[0 1]);
    82                                 md = checkfield(md,'flowequation.fe_HO','numel',[1],'values',[0 1 2]);
     82                                md = checkfield(md,'flowequation.fe_HO','numel',[1],'values',[0:3]);
    8383                                md = checkfield(md,'flowequation.fe_FS','numel',[1],'values',[0]);
    8484                                md = checkfield(md,'flowequation.borderSSA','size',[md.mesh.numberofvertices 1],'values',[0 1]);

  • issm/trunk-jpl/src/m/classes/flowequation.py

    r15621 r15636  
    6262                        md = checkfield(md,'flowequation.isFS','numel',[1],'values',[0,1])
    6363                        md = checkfield(md,'flowequation.fe_SSA','numel',[1],'values',[0,1])
    64                         md = checkfield(md,'flowequation.fe_HO','numel',[1],'values',[0,1])
     64                        md = checkfield(md,'flowequation.fe_HO','numel',[1],'values',numpy.arange(0,3+1))
    6565                        md = checkfield(md,'flowequation.fe_FS','numel',[1],'values',[0])
    6666                        md = checkfield(md,'flowequation.borderSSA','size',[md.mesh.numberofvertices],'values',[0,1])

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

    r15621 r15636  
    67016701        return StringToEnum('P2xP1')[0]
    67026702
     6703def P1xP2Enum():
     6704        """
     6705        P1XP2ENUM - Enum of P1xP2
     6706
     6707        WARNING: DO NOT MODIFY THIS FILE
     6708                                this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
     6709                                Please read src/c/shared/Enum/README for more information
     6710
     6711           Usage:
     6712              macro=P1xP2Enum()
     6713        """
     6714
     6715        return StringToEnum('P1xP2')[0]
     6716
    67036717def MINIEnum():
    67046718        """
    … …  
    79457959        """
    79467960
    7947         return 566
    7948 
     7961        return 567
     7962

  • issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m

    r15621 r15636  
    99%      macro=MaximumNumberOfEnums()
    1010
    11 macro=566;
     11macro=567;
Note: See TracChangeset for help on using the changeset viewer.