Changeset 3947


Ignore:
Timestamp:
05/26/10 09:30:24 (15 years ago)
Author:
seroussi
Message:

modified hutter solution to be done by tria or penta

Location:
issm/trunk/src/c
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/EnumDefinitions/EnumAsString.cpp

    r3828 r3947  
    161161                case NodeOnIceShelfEnum : return "NodeOnIceShelf";
    162162                case NodeOnSurfaceEnum : return "NodeOnSurface";
     163                case NumberNodeToElementConnectivityEnum : return "NumberNodeToElementConnectivity";
    163164                case PenaltyOffsetEnum : return "PenaltyOffset";
    164165                case PflagEnum : return "Pflag";

  • issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h

    r3938 r3947  
    188188        NodeOnIceShelfEnum,
    189189        NodeOnSurfaceEnum,
     190        NumberNodeToElementConnectivityEnum,
    190191        PenaltyOffsetEnum,
    191192        PflagEnum,

  • issm/trunk/src/c/EnumDefinitions/StringAsEnum.cpp

    r3828 r3947  
    159159        else if (strcmp(name,"NodeOnIceShelf")==0) return NodeOnIceShelfEnum;
    160160        else if (strcmp(name,"NodeOnSurface")==0) return NodeOnSurfaceEnum;
     161        else if (strcmp(name,"NumberNodeToElementConnectivity")==0) return NumberNodeToElementConnectivityEnum;
    161162        else if (strcmp(name,"PenaltyOffset")==0) return PenaltyOffsetEnum;
    162163        else if (strcmp(name,"Pflag")==0) return PflagEnum;

  • issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp

    r3913 r3947  
    5050        IoModelFetchData(&iomodel->rheology_B,NULL,NULL,iomodel_handle,"rheology_B");
    5151        IoModelFetchData(&iomodel->rheology_n,NULL,NULL,iomodel_handle,"rheology_n");
     52        IoModelFetchData(&iomodel->elements_type,NULL,NULL,iomodel_handle,"elements_type");
     53        IoModelFetchData(&iomodel->elementonbed,NULL,NULL,iomodel_handle,"elementonbed");
     54        IoModelFetchData(&iomodel->elementonsurface,NULL,NULL,iomodel_handle,"elementonsurface");
     55        IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater");
     56        IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
     57        IoModelFetchData(&iomodel->upperelements,NULL,NULL,iomodel_handle,"upperelements");
     58        IoModelFetchData(&iomodel->lowerelements,NULL,NULL,iomodel_handle,"lowerelements");
    5259
    5360        /*2d mesh: */
    5461        if (strcmp(iomodel->meshtype,"2d")==0){
    5562
    56                 for (i=0;i<iomodel->numberofvertices;i++){
     63                for (i=0;i<iomodel->numberofelements;i++){
    5764
    58                         if(iomodel->my_vertices[i]){
     65                        if(iomodel->my_elements[i]){
    5966
    60                                 /*Create and add penta element to elements dataset: */
    61                                 elements->AddObject(new Sing(i+1,i,iomodel));
     67                                if (*(iomodel->elements_type+2*i+0)==HutterFormulationEnum){ //create only Hutter elements
    6268
    63                                 /*Create and add material property to materials dataset: */
    64                                 materials->AddObject(new Matice(i+1,i,iomodel,1));
     69                                        /*Create and add penta element to elements dataset: */
     70                                        elements->AddObject(new Tria(i+1,i,iomodel));
    6571
     72                                        /*Create and add material property to materials dataset: */
     73                                        materials->AddObject(new Matice(i+1,i,iomodel,3));
     74                                }
    6675                        }
    67 
    6876                } //for (i=0;i<iomodel->numberofvertices;i++)
    6977        } //if (strcmp(iomodel->meshtype,"2d")==0)
    7078        else{
    7179
    72                 for (i=0;i<iomodel->numberofvertices;i++){
     80                for (i=0;i<iomodel->numberofelements;i++){
    7381
    74                         if(iomodel->my_vertices[i]){
    75                                 if(iomodel->gridonhutter[i]){
    76                                         if(!iomodel->gridonsurface[i]){
     82                        if(iomodel->my_elements[i]){
    7783
    78                                                 /*Create and add penta element to elements dataset: */
    79                                                 elements->AddObject(new Beam(i+1,i,iomodel));
     84                                if (*(iomodel->elements_type+2*i+0)==HutterFormulationEnum){ //create only Hutter elements
    8085
    81                                                 /*Create and add material property to materials dataset: */
    82                                                 materials->AddObject(new Matice(i+1,i,iomodel,2));
     86                                        /*Create and add penta element to elements dataset: */
     87                                        elements->AddObject(new Penta(i+1,i,iomodel));
    8388
    84                                         }
     89
     90                                        /*Create and add material property to materials dataset: */
     91                                        materials->AddObject(new Matice(i+1,i,iomodel,6));
     92
    8593                                }
    8694                        }
    8795
    88                 } //for (i=0;i<iomodel->numberofvertices;i++)
     96                } //for (i=0;i<iomodel->numberofelements;i++)
    8997
    9098        } //if (strcmp(iomodel->meshtype,"2d")==0)
    … …  
    92100        /*Free data: */
    93101        xfree((void**)&iomodel->elements);
     102        xfree((void**)&iomodel->elementonbed);
     103        xfree((void**)&iomodel->elementonsurface);
     104        xfree((void**)&iomodel->elements_type);
    94105        xfree((void**)&iomodel->gridonhutter);
    95106        xfree((void**)&iomodel->thickness);
    96107        xfree((void**)&iomodel->surface);
    97108        xfree((void**)&iomodel->bed);
    98         xfree((void**)&iomodel->gridonsurface);
     109        xfree((void**)&iomodel->gridonbed);
     110        xfree((void**)&iomodel->elementonwater);
    99111        xfree((void**)&iomodel->uppernodes);
    100112        xfree((void**)&iomodel->drag_coefficient);
    101113        xfree((void**)&iomodel->rheology_B);
    102114        xfree((void**)&iomodel->rheology_n);
     115        xfree((void**)&iomodel->upperelements);
     116        xfree((void**)&iomodel->lowerelements);
    103117
    104118        /*Add new constrant material property to materials, at the end: */
    105119        if (strcmp(iomodel->meshtype,"2d")==0){
    106                 materials->AddObject(new Matpar(iomodel->numberofvertices+1,iomodel));                          //put it at the end of the materials
     120                materials->AddObject(new Matpar(iomodel->numberofelements+1,iomodel));                          //put it at the end of the materials
    107121        }
    108122        else{
    109                 materials->AddObject(new Matpar(iomodel->numberofvertices2d*(iomodel->numlayers-1)+1,iomodel)); //put it at the end of the materials
     123                materials->AddObject(new Matpar(iomodel->numberofelements+1,iomodel)); //put it at the end of the materials
    110124        }
    111125               
    … …  
    125139        IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet");
    126140        IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf");
     141        IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
     142        CreateNumberNodeToElementConnectivity(iomodel);
    127143
    128144        for (i=0;i<iomodel->numberofvertices;i++){
    … …  
    153169        xfree((void**)&iomodel->gridonicesheet);
    154170        xfree((void**)&iomodel->gridoniceshelf);
     171        xfree((void**)&iomodel->elements);
     172        xfree((void**)&iomodel->numbernodetoelementconnectivity);
    155173
    156174        /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these

  • issm/trunk/src/c/objects/Elements/Beam.cpp

    r3944 r3947  
    414414        int       numberofdofspernode;
    415415        bool onbed;
    416        
    417        
     416        bool onsurface;
     417        int connectivity[2];
     418        double one0,one1;
     419       
     420        /*dynamic objects pointed to by hooks: */
     421        Node**  nodes=NULL;
     422
     423        /*recover objects from hooks: */
     424        nodes=(Node**)hnodes.deliverp();
     425
     426        connectivity[0]=nodes[0]->GetConnectivity();
     427        connectivity[1]=nodes[1]->GetConnectivity();
     428       
     429        one0=1/(double)connectivity[0];
     430        one1=1/(double)connectivity[1];
    418431        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     432        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
    419433
    420434        GetDofList(&doflist[0],&numberofdofspernode);
    421435
    422436        if (onbed){
    423                 Ke_gg[0][0]=1;
    424                 Ke_gg[1][1]=1;
    425                 Ke_gg[2][0]=-1;
    426                 Ke_gg[2][2]=1;
    427                 Ke_gg[3][1]=-1;
    428                 Ke_gg[3][3]=1;
    429         }
    430         else{
    431                 Ke_gg[2][0]=-1;
    432                 Ke_gg[2][2]=1;
    433                 Ke_gg[3][1]=-1;
    434                 Ke_gg[3][3]=1;
     437                Ke_gg[0][0]=one0;
     438                Ke_gg[1][1]=one0;
     439                Ke_gg[2][0]=-2*one1;
     440                Ke_gg[2][2]=2*one1;
     441                Ke_gg[3][1]=-2*one1;
     442                Ke_gg[3][3]=2*one1;
     443        }
     444        else if (onsurface){
     445                Ke_gg[2][0]=-one1;
     446                Ke_gg[2][2]=one1;
     447                Ke_gg[3][1]=-one1;
     448                Ke_gg[3][3]=one1;
     449        }
     450        else{ //node is on two horizontal layers and beams include the values only once, so the have to use half of the connectivity
     451                Ke_gg[2][0]=-2*one1;
     452                Ke_gg[2][2]=2*one1;
     453                Ke_gg[3][1]=-2*one1;
     454                Ke_gg[3][3]=2*one1;
    435455        }
    436456
    … …  
    495515       
    496516        bool onbed;
     517        bool onsurface;
     518        int  connectivity[2];
    497519
    498520        /*dynamic objects pointed to by hooks: */
    … …  
    511533        /*recover some inputs: */
    512534        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     535        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
    513536
    514537        /*recover parameters: */
    … …  
    524547        inputs->GetParameterValue(&thickness,&gauss1[0],ThicknessEnum);
    525548
     549        connectivity[0]=nodes[0]->GetConnectivity();
     550        connectivity[1]=nodes[1]->GetConnectivity();
     551
    526552        //Get all element grid data:
    527553        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    … …  
    551577                GetJacobianDeterminant(&Jdet, &z_list[0],gauss_coord);
    552578               
    553                 for(j=0;j<NDOF2;j++){
    554                         pe_g_gaussian[NDOF2+j]=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss_weight;
     579                /*Add contribution*/
     580                if (onsurface){
     581                        for(j=0;j<NDOF2;j++){
     582                                pe_g_gaussian[NDOF2+j]=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss_weight/(double)connectivity[1];
     583                        }
     584                }
     585                else{//connectivity is too large, should take only half on it
     586                        for(j=0;j<NDOF2;j++){
     587                                pe_g_gaussian[NDOF2+j]=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss_weight*2/(double)connectivity[1];
     588                        }
    555589                }
    556590               
    … …  
    570604
    571605                //Add to pe:
    572                 pe_g[0]+=ub;
    573                 pe_g[1]+=vb;
     606                pe_g[0]+=ub/(double)connectivity[0];
     607                pe_g[1]+=vb/(double)connectivity[0];
    574608        }
    575609

  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r3946 r3947  
    5353
    5454        /*hooks: */
     55        ISSMASSERT(iomodel->elements);
    5556        for(i=0;i<6;i++){ //go recover node ids, needed to initialize the node hook.
    5657                penta_node_ids[i]=(int)*(iomodel->elements+6*index+i); //ids for vertices are in the elements array from Matlab
    … …  
    5960        penta_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object
    6061
     62        ISSMASSERT(iomodel->upperelements);
    6163        if isnan(iomodel->upperelements[index]){
    6264                penta_elements_ids[0]=this->id; //upper penta is the same penta
    … …  
    6567                penta_elements_ids[0]=(int)(iomodel->upperelements[index]);
    6668        }
     69        ISSMASSERT(iomodel->lowerelements);
    6770        if isnan(iomodel->lowerelements[index]){
    6871                penta_elements_ids[1]=this->id; //lower penta is the same penta
    … …  
    7174                penta_elements_ids[1]=(int)(iomodel->lowerelements[index]);
    7275        }
    73 
    7476        this->InitHookNodes(penta_node_ids); this->nodes=NULL;
    7577        this->InitHookMatice(penta_matice_id);this->matice=NULL;
    … …  
    7981        //intialize inputs, and add as many inputs per element as requested:
    8082        this->inputs=new Inputs();
    81 
    8283        if (iomodel->thickness) {
    8384                for(i=0;i<6;i++)nodeinputs[i]=iomodel->thickness[penta_node_ids[i]-1];
    … …  
    125126                this->inputs->AddInput(new PentaVertexInput(DhDtEnum,nodeinputs));
    126127        }
    127 
    128128        /*vx,vy and vz: */
    129129        if (iomodel->vx) {
    … …  
    178178                this->inputs->AddInput(new PentaVertexInput(VzOldEnum,nodeinputs));
    179179        }
    180 
    181180        if (iomodel->elementoniceshelf) this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,(IssmBool)iomodel->elementoniceshelf[index]));
    182181        if (iomodel->elementonbed) this->inputs->AddInput(new BoolInput(ElementOnBedEnum,(IssmBool)iomodel->elementonbed[index]));
    … …  
    250249        /*point parameters to real dataset: */
    251250        this->parameters=parametersin;
    252 
    253251}
    254252/*}}}*/
    … …  
    330328                                name==SurfaceEnum ||
    331329                                name==BedEnum ||
     330                                name==SurfaceSlopexEnum ||
     331                                name==SurfaceSlopeyEnum ||
    332332                                name==MeltingRateEnum ||
    333333                                name==AccumulationRateEnum ||

  • issm/trunk/src/c/objects/Elements/Sing.cpp

    r3944 r3947  
    333333        int    doflist[numdofs];
    334334        int    numberofdofspernode;
     335        int    connectivity;
     336
     337        /*dynamic objects pointed to by hooks: */
     338        Node**  nodes=NULL;
     339
     340        /*recover objects from hooks: */
     341        nodes=(Node**)hnodes.delivers();
     342
     343        /*Find connectivity of the node and divide Ke_gg by this connectivity*/
     344        connectivity=nodes[0]->GetConnectivity();
     345        Ke_gg[0][0]=1/(double)connectivity;
     346        Ke_gg[1][1]=1/(double)connectivity;
    335347
    336348        GetDofList(&doflist[0],&numberofdofspernode);
    … …  
    374386        double    rho_ice,gravity,n,B;
    375387        double    thickness;
     388        int       connectivity;
    376389
    377390        /*dynamic objects pointed to by hooks: */
    378         Node**  node=NULL;
     391        Node**  nodes=NULL;
    379392        Matpar* matpar=NULL;
    380393        Matice* matice=NULL;
    381394
    382395        /*recover objects from hooks: */
    383         node=(Node**)hnodes.deliverp();
     396        nodes=(Node**)hnodes.deliverp();
    384397        matpar=(Matpar*)hmatpar.delivers();
    385398        matice=(Matice*)hmatice.delivers();
    … …  
    389402
    390403        GetDofList(&doflist[0],&numberofdofspernode);
     404
     405        //Get connectivity of the node
     406        connectivity=nodes[0]->GetConnectivity();
    391407
    392408        //compute slope2
    … …  
    406422        constant_part=-2*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2));
    407423
    408         pe_g[0]=ub-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[0];
    409         pe_g[1]=vb-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[1];
     424        pe_g[0]=(ub-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[0])/(double)connectivity;
     425        pe_g[1]=(vb-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[1])/(double)connectivity;
    410426
    411427        VecSetValues(pg,numdofs,doflist,(const double*)pe_g,ADD_VALUES);
    412 
    413428
    414429}

  • issm/trunk/src/c/objects/Node.cpp

    r3858 r3947  
    165165        if (iomodel->gridoniceshelf) this->inputs->AddInput(new BoolInput(NodeOnIceShelfEnum,(IssmBool)iomodel->gridoniceshelf[i]));
    166166        if (iomodel->gridonicesheet) this->inputs->AddInput(new BoolInput(NodeOnIceSheetEnum,(IssmBool)iomodel->gridonicesheet[i]));
     167        NumberNodeToElementConnectivityEnum;
     168        if (iomodel->numbernodetoelementconnectivity) this->inputs->AddInput(new IntInput(NumberNodeToElementConnectivityEnum,iomodel->numbernodetoelementconnectivity[i]));
    167169
    168170}
    … …  
    771773}
    772774/*}}}*/
     775/*FUNCTION Node::GetConnectivity {{{2*/
     776int Node::GetConnectivity(){
     777        int connectivity;
     778
     779        /*recover parameters: */
     780        inputs->GetParameterValue(&connectivity,NumberNodeToElementConnectivityEnum);
     781
     782        return connectivity;
     783}
     784/*}}}*/
    773785/*FUNCTION Node::GetNumberOfDofs{{{2*/
    774786int   Node::GetNumberOfDofs(){

  • issm/trunk/src/c/objects/Node.h

    r3858 r3947  
    7070                int   GetDof(int dofindex);
    7171                void  CreateVecSets(Vec pv_g,Vec pv_m,Vec pv_n,Vec pv_f,Vec pv_s);
     72                int   GetConnectivity();
    7273                void  GetDofList(int* outdoflist,int* pnumberofdofspernode);
    7374                int   GetDofList1(void);
Note: See TracChangeset for help on using the changeset viewer.