Changeset 18173


Ignore:
Timestamp:
06/23/14 09:47:02 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: added CrouzeixRaviart elements in FS

Location:
issm/trunk-jpl/src
Files:
11 edited

Legend:

Unmodified
Added
Removed

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

    r18151 r18173  
    414414                        /*Deduce velocity interpolation from finite element*/
    415415                        switch(finiteelement){
    416                                 case P1P1Enum          : finiteelement = P1Enum;       break;
    417                                 case P1P1GLSEnum       : finiteelement = P1Enum;       break;
    418                                 case MINIcondensedEnum : finiteelement = P1bubbleEnum; break;
    419                                 case MINIEnum          : finiteelement = P1bubbleEnum; break;
    420                                 case TaylorHoodEnum    : finiteelement = P2Enum;       break;
    421                                 case XTaylorHoodEnum   : finiteelement = P2Enum;       break;
    422                                 case OneLayerP4zEnum   : finiteelement = P2xP4Enum;    break;
     416                                case P1P1Enum            : finiteelement = P1Enum;       break;
     417                                case P1P1GLSEnum         : finiteelement = P1Enum;       break;
     418                                case MINIcondensedEnum   : finiteelement = P1bubbleEnum; break;
     419                                case MINIEnum            : finiteelement = P1bubbleEnum; break;
     420                                case TaylorHoodEnum      : finiteelement = P2Enum;       break;
     421                                case XTaylorHoodEnum     : finiteelement = P2Enum;       break;
     422                                case OneLayerP4zEnum     : finiteelement = P2xP4Enum;    break;
     423                                case CrouzeixRaviartEnum : finiteelement = P2bubbleEnum; break;
    423424                                default: _error_("finite element "<<EnumToStringx(finiteelement)<<" not supported");
    424425                        }
    … …  
    494495                                                        if(reCast<int,IssmDouble>(vertices_type[i])==NoneApproximationEnum){
    495496                                                                constraints->AddObject(new SpcStatic(count+1,iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+i+1,0,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum));
     497                                                                count++;
     498                                                        }
     499                                                }
     500                                        }
     501                                        break;
     502                                case P2bubbleEnum:
     503                                        for(i=0;i<iomodel->numberofvertices;i++){
     504                                                if(iomodel->my_vertices[i]){
     505                                                        if(reCast<int,IssmDouble>(vertices_type[i])==NoneApproximationEnum){
     506                                                                constraints->AddObject(new SpcStatic(count+1,iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements+iomodel->numberoffaces+i+1,0,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum));
    496507                                                                count++;
    497508                                                        }

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

    r18171 r18173  
    23622362                        penta_node_ids[35]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->numberoffaces+iomodel->elements[6*index+5];
    23632363                        break;
     2364                case CrouzeixRaviartEnum:
     2365                        numnodes         = 25;
     2366                        penta_node_ids   = xNew<int>(numnodes);
     2367                        penta_node_ids[ 0]=iomodel->nodecounter+iomodel->elements[6*index+0];
     2368                        penta_node_ids[ 1]=iomodel->nodecounter+iomodel->elements[6*index+1];
     2369                        penta_node_ids[ 2]=iomodel->nodecounter+iomodel->elements[6*index+2];
     2370                        penta_node_ids[ 3]=iomodel->nodecounter+iomodel->elements[6*index+3];
     2371                        penta_node_ids[ 4]=iomodel->nodecounter+iomodel->elements[6*index+4];
     2372                        penta_node_ids[ 5]=iomodel->nodecounter+iomodel->elements[6*index+5];
     2373                        penta_node_ids[ 6]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+0]+1;
     2374                        penta_node_ids[ 7]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+1]+1;
     2375                        penta_node_ids[ 8]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+2]+1;
     2376                        penta_node_ids[ 9]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+3]+1;
     2377                        penta_node_ids[10]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+4]+1;
     2378                        penta_node_ids[11]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+5]+1;
     2379                        penta_node_ids[12]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+6]+1;
     2380                        penta_node_ids[13]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1;
     2381                        penta_node_ids[14]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1;
     2382                        penta_node_ids[15]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+2]+1;
     2383                        penta_node_ids[16]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+3]+1;
     2384                        penta_node_ids[17]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+4]+1;
     2385                        penta_node_ids[18]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+index+1;
     2386
     2387                        penta_node_ids[19]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->numberofelements+6*index+1;
     2388                        penta_node_ids[20]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->numberofelements+6*index+2;
     2389                        penta_node_ids[21]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->numberofelements+6*index+3;
     2390                        penta_node_ids[22]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->numberofelements+6*index+4;
     2391                        penta_node_ids[23]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->numberofelements+6*index+5;
     2392                        penta_node_ids[24]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->numberofelements+6*index+6;
     2393                        break;
    23642394                default:
    23652395                        _error_("Finite element "<<EnumToStringx(finiteelement_type)<<" not supported yet");

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

    r18171 r18173  
    947947                case P0Enum:                return NUMNODESP0;
    948948                case P1Enum:                return NUMNODESP1;
     949                case P1DGEnum:              return NUMNODESP1;
    949950                case P1bubbleEnum:          return NUMNODESP1b;
    950951                case P1bubblecondensedEnum: return NUMNODESP1b;
    … …  
    954955                case P2xP1Enum:             return NUMNODESP2xP1;
    955956                case P1xP2Enum:             return NUMNODESP1xP2;
     957                case P2xP4Enum:             return NUMNODESP2xP4;
     958                case P1xP3Enum:             return NUMNODESP1xP3;
    956959                case P1P1Enum:              return NUMNODESP1*2;
    957960                case P1P1GLSEnum:           return NUMNODESP1*2;
    … …  
    960963                case TaylorHoodEnum:        return NUMNODESP2+NUMNODESP1;
    961964                case OneLayerP4zEnum:       return NUMNODESP2xP4+NUMNODESP1;
    962                 case P2xP4Enum:             return NUMNODESP2xP4;
    963                 case P1xP3Enum:             return NUMNODESP1xP3;
     965                case CrouzeixRaviartEnum:   return NUMNODESP2b+NUMNODESP1;
    964966                default:       _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
    965967        }
    … …  
    971973
    972974        switch(fe_stokes){
    973                 case P1P1Enum:          return P1Enum;
    974                 case P1P1GLSEnum:       return P1Enum;
    975                 case MINIcondensedEnum: return P1bubbleEnum;
    976                 case MINIEnum:          return P1bubbleEnum;
    977                 case TaylorHoodEnum:    return P2Enum;
    978                 case OneLayerP4zEnum:   return P2xP4Enum;
     975                case P1P1Enum:           return P1Enum;
     976                case P1P1GLSEnum:        return P1Enum;
     977                case MINIcondensedEnum:  return P1bubbleEnum;
     978                case MINIEnum:           return P1bubbleEnum;
     979                case TaylorHoodEnum:     return P2Enum;
     980                case OneLayerP4zEnum:    return P2xP4Enum;
     981                case CrouzeixRaviartEnum:return P2bubbleEnum;
    979982                default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
    980983        }
    … …  
    986989
    987990        switch(fe_stokes){
    988                 case P1P1Enum:          return P1Enum;
    989                 case P1P1GLSEnum:       return P1Enum;
    990                 case MINIcondensedEnum: return P1Enum;
    991                 case MINIEnum:          return P1Enum;
    992                 case TaylorHoodEnum:    return P1Enum;
    993                 case OneLayerP4zEnum:   return P1Enum;
     991                case P1P1Enum:           return P1Enum;
     992                case P1P1GLSEnum:        return P1Enum;
     993                case MINIcondensedEnum:  return P1Enum;
     994                case MINIEnum:           return P1Enum;
     995                case TaylorHoodEnum:     return P1Enum;
     996                case OneLayerP4zEnum:    return P1Enum;
     997                case CrouzeixRaviartEnum:return P1DGEnum;
    994998                default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
    995999        }
    … …  
    10541058                        break;
    10551059                case P2Enum:
     1060                        numindices = 6;
     1061                        indices    = xNew<int>(numindices);
     1062                        indices[0] = 0;
     1063                        indices[1] = 1;
     1064                        indices[2] = 2;
     1065                        indices[3] = 9;
     1066                        indices[4] = 10;
     1067                        indices[5] = 11;
     1068                        break;
     1069                case P2bubbleEnum:
    10561070                        numindices = 6;
    10571071                        indices    = xNew<int>(numindices);

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

    r18170 r18173  
    19751975                        tria_node_ids[8]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[3*index+2];
    19761976                        break;
     1977                case CrouzeixRaviartEnum:
     1978                        numnodes        = 10;
     1979                        tria_node_ids   = xNew<int>(numnodes);
     1980                        tria_node_ids[0]=iomodel->nodecounter+iomodel->elements[3*index+0];
     1981                        tria_node_ids[1]=iomodel->nodecounter+iomodel->elements[3*index+1];
     1982                        tria_node_ids[2]=iomodel->nodecounter+iomodel->elements[3*index+2];
     1983                        tria_node_ids[3]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+0]+1;
     1984                        tria_node_ids[4]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+1]+1;
     1985                        tria_node_ids[5]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+2]+1;
     1986                        tria_node_ids[6]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+index+1;
     1987
     1988                        tria_node_ids[7]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements+3*index+1;
     1989                        tria_node_ids[8]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements+3*index+2;
     1990                        tria_node_ids[9]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements+3*index+3;
     1991                        break;
    19771992                default:
    19781993                        _error_("Finite element "<<EnumToStringx(finiteelement_type)<<" not supported yet");

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

    r18170 r18173  
    420420                case TaylorHoodEnum:        return NUMNODESP2+NUMNODESP1;
    421421                case XTaylorHoodEnum:       return NUMNODESP2+NUMNODESP1;
     422                case CrouzeixRaviartEnum:   return NUMNODESP2b+NUMNODESP1;
    422423                default: _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
    423424        }
    … …  
    429430
    430431        switch(fe_stokes){
    431                 case P1P1Enum:          return P1Enum;
    432                 case P1P1GLSEnum:       return P1Enum;
    433                 case MINIcondensedEnum: return P1bubbleEnum;
    434                 case MINIEnum:          return P1bubbleEnum;
    435                 case TaylorHoodEnum:    return P2Enum;
    436                 case XTaylorHoodEnum:   return P2Enum;
     432                case P1P1Enum:           return P1Enum;
     433                case P1P1GLSEnum:        return P1Enum;
     434                case MINIcondensedEnum:  return P1bubbleEnum;
     435                case MINIEnum:           return P1bubbleEnum;
     436                case TaylorHoodEnum:     return P2Enum;
     437                case XTaylorHoodEnum:    return P2Enum;
     438                case CrouzeixRaviartEnum:return P2bubbleEnum;
    437439                default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
    438440        }
    … …  
    444446
    445447        switch(fe_stokes){
    446                 case P1P1Enum:          return P1Enum;
    447                 case P1P1GLSEnum:       return P1Enum;
    448                 case MINIcondensedEnum: return P1Enum;
    449                 case MINIEnum:          return P1Enum;
    450                 case TaylorHoodEnum:    return P1Enum;
    451                 case XTaylorHoodEnum:   return P1Enum;
     448                case P1P1Enum:            return P1Enum;
     449                case P1P1GLSEnum:         return P1Enum;
     450                case MINIcondensedEnum:   return P1Enum;
     451                case MINIEnum:            return P1Enum;
     452                case TaylorHoodEnum:      return P1Enum;
     453                case XTaylorHoodEnum:     return P1Enum;
     454                case CrouzeixRaviartEnum: return P1DGEnum;
    452455                default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
    453456        }

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

    r18170 r18173  
    1313        /*Intermediaries*/
    1414        int   i,j,counter,vnodes,lid=0;
    15         int   numberoffaces;
     15        int   numberoffaces,elementnbv;
    1616        int   id0 = iomodel->nodecounter;
    1717        bool *my_faces = NULL;
    … …  
    442442                        }
    443443                        break;
     444                case CrouzeixRaviartEnum:
     445                        _assert_(approximation==FSApproximationEnum);
     446                        /*P2b velocity*/
     447                        EdgesPartitioning(&my_edges,iomodel);
     448                        for(i=0;i<iomodel->numberofvertices;i++){
     449                                if(iomodel->my_vertices[i]){
     450                                        nodes->AddObject(new Node(id0+i+1,i,lid++,i,iomodel,analysis,FSvelocityEnum));
     451                                }
     452                        }
     453                        for(i=0;i<iomodel->numberofedges;i++){
     454                                if(my_edges[i]){
     455                                        nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,iomodel,analysis,FSvelocityEnum));
     456                                }
     457                        }
     458                        id0 = id0+iomodel->numberofvertices+iomodel->numberofedges;
     459                        if(iomodel->meshelementtype==PentaEnum){
     460                                FacesPartitioning(&my_faces,iomodel);
     461                                for(i=0;i<iomodel->numberoffaces;i++){
     462                                        if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
     463                                                if(my_faces[i]){
     464                                                        node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,iomodel,analysis,FSvelocityEnum);
     465                                                        nodes->AddObject(node);
     466                                                }
     467                                        }
     468                                        else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/
     469                                                /*Nothing*/
     470                                        }
     471                                        else{
     472                                                _error_("not supported");
     473                                        }
     474                                }
     475                                id0 = id0+iomodel->numberoffaces;
     476                        }
     477                        for(i=0;i<iomodel->numberofelements;i++){
     478                                if(iomodel->my_elements[i]){
     479                                        nodes->AddObject(new Node(id0+i+1,id0-iomodel->nodecounter+i,lid++,0,iomodel,analysis,FSvelocityEnum));
     480                                }
     481                        }
     482
     483                        /*P1 DG pressure*/
     484                        vnodes = id0+iomodel->numberofelements;
     485                        switch(iomodel->meshelementtype){
     486                                case TriaEnum:  elementnbv = 3; break;
     487                                case TetraEnum: elementnbv = 4; break;
     488                                case PentaEnum: elementnbv = 6; break;
     489                                default:        _error_("mesh dimension not supported yet");
     490                        }
     491                        for(i=0;i<iomodel->numberofelements;i++){
     492                                if(iomodel->my_elements[i]){
     493                                        for(j=0;j<elementnbv;j++){
     494                                                nodes->AddObject(new Node(vnodes+elementnbv*i+j+1,vnodes-iomodel->nodecounter+elementnbv*i+j,lid++,iomodel->elements[+elementnbv*i+j]-1,iomodel,analysis,FSpressureEnum));
     495
     496                                        }
     497                                }
     498                        }
     499                        break;
    444500
    445501                default:

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

    r18170 r18173  
    615615        XTaylorHoodEnum,
    616616        OneLayerP4zEnum,
     617        CrouzeixRaviartEnum,
    617618        /*}}}*/
    618619        /*Results{{{*/

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

    r18170 r18173  
    604604                case XTaylorHoodEnum : return "XTaylorHood";
    605605                case OneLayerP4zEnum : return "OneLayerP4z";
     606                case CrouzeixRaviartEnum : return "CrouzeixRaviart";
    606607                case SaveResultsEnum : return "SaveResults";
    607608                case BoolExternalResultEnum : return "BoolExternalResult";

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

    r18170 r18173  
    616616              else if (strcmp(name,"XTaylorHood")==0) return XTaylorHoodEnum;
    617617              else if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum;
     618              else if (strcmp(name,"CrouzeixRaviart")==0) return CrouzeixRaviartEnum;
    618619              else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum;
    619620              else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
    … …  
    628629              else if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum;
    629630              else if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum;
    630               else if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"MassfluxatgateName")==0) return MassfluxatgateNameEnum;
     634              if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum;
     635              else if (strcmp(name,"MassfluxatgateName")==0) return MassfluxatgateNameEnum;
    635636              else if (strcmp(name,"MassfluxatgateSegments")==0) return MassfluxatgateSegmentsEnum;
    636637              else if (strcmp(name,"MisfitName")==0) return MisfitNameEnum;

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

    r18171 r18173  
    134134                                md = checkfield(md,'fieldname','flowequation.fe_SSA','values',{'P1','P1bubble','P1bubblecondensed','P2','P2bubble'});
    135135                                md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',{'P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2','P2bubble','P1xP3','P2xP4'});
    136                                 md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',{'P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','XTaylorHood','OneLayerP4z'});
     136                                md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',{'P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','XTaylorHood','OneLayerP4z','CrouzeixRaviart'});
    137137                                md = checkfield(md,'fieldname','flowequation.XTH_r','numel',[1],'>',0.);
    138138                                md = checkfield(md,'fieldname','flowequation.XTH_theta','numel',[1],'>=',0.,'<',0.5);

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

    r18171 r18173  
    7878                        md = checkfield(md,'fieldname','flowequation.fe_SSA','values',['P1','P1bubble','P1bubblecondensed','P2','P2bubble'])
    7979                        md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',['P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2','P2bubble','P1xP3','P2xP4'])
    80                         md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',['P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','XTaylorHood','OneLayerP4z'])
     80                        md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',['P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','XTaylorHood','OneLayerP4z','CrouzeixRaviart'])
    8181                        md = checkfield(md,'fieldname','flowequation.borderSSA','size',[md.mesh.numberofvertices],'values',[0,1])
    8282                        md = checkfield(md,'fieldname','flowequation.borderHO','size',[md.mesh.numberofvertices],'values',[0,1])
Note: See TracChangeset for help on using the changeset viewer.