Changeset 5935


Ignore:
Timestamp:
09/22/10 07:47:21 (15 years ago)
Author:
Mathieu Morlighem
Message:

more ElementVector

Location:
issm/trunk/src/c/objects/Loads
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/objects/Loads/Icefront.cpp

    r5911 r5935  
    325325
    326326        /*Retrieve parameters: */
     327        ElementVector* pe=NULL;
    327328        int analysis_type;
    328329        this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    329330
    330331        /*Just branch to the correct element icefront vector generator, according to the type of analysis we are carrying out: */
    331         if (analysis_type==DiagnosticHorizAnalysisEnum){
    332                 int type;
    333                 inputs->GetParameterValue(&type,TypeEnum);
    334                 if (type==MacAyeal2dIceFrontEnum){
    335                         CreatePVectorDiagnosticMacAyeal2d( pg,pf);
    336                 }
    337                 else if (type==MacAyeal3dIceFrontEnum){
    338                         CreatePVectorDiagnosticMacAyeal3d( pg);
    339                 }
    340                 else if (type==PattynIceFrontEnum){
    341                         CreatePVectorDiagnosticPattyn( pg);
    342                 }
    343                 else if (type==StokesIceFrontEnum){
    344                         CreatePVectorDiagnosticStokes( pg);
    345                 }
    346                 else ISSMERROR("Icefront type %s not supported yet",EnumToString(type));
    347         }
    348         else if (analysis_type==AdjointHorizAnalysisEnum){
    349                 return;
    350         }
    351         else{
    352                 ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type));
     332        switch(analysis_type){
     333                case DiagnosticHorizAnalysisEnum:
     334                        pe=CreatePVectorDiagnosticHoriz();
     335                        break;
     336                case AdjointHorizAnalysisEnum:
     337                        pe=CreatePVectorAdjointHoriz();
     338                        break;
     339                default:
     340                        ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type));
     341        }
     342
     343        /*Add to global Vector*/
     344        if(pe){
     345                pe->AddToGlobal(pg,pf);
     346                delete pe;
    353347        }
    354348}
    … …  
    424418
    425419/*Icefront numerics: */
     420/*FUNCTION Icefront::CreatePVectorDiagnosticHoriz {{{1*/
     421ElementVector* Icefront::CreatePVectorDiagnosticHoriz(void){
     422
     423        int type;
     424        inputs->GetParameterValue(&type,TypeEnum);
     425
     426        switch(type){
     427                case MacAyeal2dIceFrontEnum:
     428                        return CreatePVectorDiagnosticMacAyeal2d();
     429                case MacAyeal3dIceFrontEnum:
     430                        return CreatePVectorDiagnosticMacAyeal3d();
     431                case PattynIceFrontEnum:
     432                        return CreatePVectorDiagnosticPattyn();
     433                case StokesIceFrontEnum:
     434                        return CreatePVectorDiagnosticStokes();
     435                default:
     436                        ISSMERROR("Icefront type %s not supported yet",EnumToString(type));
     437        }
     438}
     439/*}}}*/
     440/*FUNCTION Icefront::CreatePVectorAdjointHoriz {{{1*/
     441ElementVector* Icefront::CreatePVectorAdjointHoriz(void){
     442
     443        /*No load vector applied to the adjoint*/
     444        return NULL;
     445}
     446/*}}}*/
    426447/*FUNCTION Icefront::CreatePVectorDiagnosticMacAyeal2d{{{1*/
    427 void Icefront::CreatePVectorDiagnosticMacAyeal2d( Vec pg,Vec pf){
     448ElementVector* Icefront::CreatePVectorDiagnosticMacAyeal2d(void){
    428449
    429450        /*Constants*/
    430451        const int numnodes= NUMVERTICESSEG;
    431452        const int numdofs = numnodes *NDOF2;
    432 
    433         /*output: */
    434         ElementVector* pe=NULL;
    435453
    436454        /*Intermediary*/
    … …  
    440458        double     water_pressure,air_pressure,surface_under_water,base_under_water;
    441459        double     xyz_list[numnodes][3];
    442         double     pe_g[numdofs] = {0.0};
    443460        double     normal[2];
    444461        double     L[2];
    … …  
    447464        Tria* tria=((Tria*)element);
    448465
    449         /* Get dof list and node coordinates: */
     466        /*Initialize Element vector and return if necessary*/
     467        if(tria->IsOnWater()) return NULL;
     468        ElementVector* pe=NewElementVector(nodes,NUMVERTICESSEG,this->parameters);
     469
     470        /*Retrieve all inputs and parameters*/
    450471        GetVerticesCoordinates(&xyz_list[0][0],nodes,numnodes);
    451 
    452         /*Recover inputs and parameters */
    453472        Input* thickness_input=tria->inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input);
    454473        Input* bed_input      =tria->inputs->GetInput(BedEnum);       ISSMASSERT(bed_input);
    455474        inputs->GetParameterValue(&fill,FillEnum);
    456 
    457475        rho_water=matpar->GetRhoWater();
    458476        rho_ice  =matpar->GetRhoIce();
    459477        gravity  =matpar->GetG();
    460 
    461478        GetSegmentNormal(&normal[0],xyz_list);
    462479
     480        /*Start looping on Gaussian points*/
    463481        index1=tria->GetNodeIndex(nodes[0]);
    464482        index2=tria->GetNodeIndex(nodes[1]);
    … …  
    492510
    493511                for (int i=0;i<numnodes;i++){
    494                         pe_g[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*L[i];
    495                         pe_g[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*L[i];
     512                        pe->values[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*L[i];
     513                        pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*L[i];
    496514                }
    497515        }
    498516
    499         /*Initialize element matrix: */
    500         pe=NewElementVector(nodes,NUMVERTICESSEG,this->parameters,MacAyealApproximationEnum);
    501 
    502         /*Add pe_g values to pe element stifness load: */
    503         pe->AddValues(&pe_g[0]);
    504 
    505         /*Add pe element load vector to global load vector: */
    506         pe->AddToGlobal(pg,pf);
    507 
    508         /*Free ressources:*/
    509         delete pe;
    510 
    511         /*Free ressources:*/
     517        /*Clean up and return*/
    512518        delete gauss;
     519        return pe;
    513520}
    514521/*}}}*/
    515522/*FUNCTION Icefront::CreatePVectorDiagnosticMacAyeal3d{{{1*/
    516 void Icefront::CreatePVectorDiagnosticMacAyeal3d( Vec pg){
     523ElementVector* Icefront::CreatePVectorDiagnosticMacAyeal3d(void){
    517524
    518525        Icefront* icefront=NULL;
    … …  
    525532
    526533        /*Return if not on bed*/
    527         penta->inputs->GetParameterValue(&onbed,ElementOnBedEnum);
    528         if(!onbed) return;
     534        if(!penta->IsOnBed() || penta->IsOnWater()) return NULL;
    529535
    530536        /*Spawn Tria and call MacAyeal2d*/
    … …  
    533539        icefront->element=tria;
    534540        icefront->inputs->AddInput(new IntInput(TypeEnum,MacAyeal2dIceFrontEnum));
    535         icefront->CreatePVectorDiagnosticMacAyeal2d(pg,NULL);
    536 
    537         /*clean-up*/
     541        ElementVector* pe=icefront->CreatePVectorDiagnosticMacAyeal2d();
     542
     543        /*clean-up and return*/
    538544        delete tria->matice;
    539545        delete tria;
    540546        delete icefront;
     547        return pe;
    541548}
    542549/*}}}*/
    543550/*FUNCTION Icefront::CreatePVectorDiagnosticPattyn{{{1*/
    544 void Icefront::CreatePVectorDiagnosticPattyn( Vec pg){
     551ElementVector* Icefront::CreatePVectorDiagnosticPattyn(void){
    545552
    546553        /*Constants*/
    … …  
    555562        double      xyz_list[NUMVERTICESQUA][3];
    556563        double      normal[3];
    557         double      pe_g[numdofs] = {0.0};
    558564        double      l1l4[4];
    559         int        *doflist = NULL;
    560565        GaussPenta *gauss = NULL;
    561         Penta      *penta = NULL;
    562 
    563         /* Get dof list and node coordinates: */
    564         penta=(Penta*)element;
    565         GetDofList(&doflist,PattynApproximationEnum,GsetEnum);
     566
     567        Penta* penta=(Penta*)element;
     568
     569        /*Initialize Element vector and return if necessary*/
     570        if(penta->IsOnWater()) return NULL;
     571        ElementVector* pe=NewElementVector(nodes,NUMVERTICESQUA,this->parameters);
     572
     573        /*Retrieve all inputs and parameters*/
    566574        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICESQUA);
    567 
    568         /*Recover inputs and parameters */
    569575        Input* surface_input  =penta->inputs->GetInput(SurfaceEnum);   ISSMASSERT(surface_input);
    570576        inputs->GetParameterValue(&fill,FillEnum);
    … …  
    608614                pressure = ice_pressure + water_pressure + air_pressure;
    609615
    610                 for(i=0;i<NUMVERTICESQUA;i++){
    611                         for(j=0;j<NDOF2;j++){
    612                                 pe_g[i*NDOF2+j]+=Jdet*gauss->weight*pressure*l1l4[i]*normal[j];
    613                         }
    614                 }
    615         }
    616 
    617         VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES);
    618 
    619         /*Free ressources:*/
     616                for(i=0;i<NUMVERTICESQUA;i++) for(j=0;j<NDOF2;j++) pe->values[i*NDOF2+j]+=Jdet*gauss->weight*pressure*l1l4[i]*normal[j];
     617        }
     618
     619        /*Clean up and return*/
    620620        delete gauss;
    621         xfree((void**)&doflist);
     621        return pe;
    622622}
    623623/*}}}*/
    624624/*FUNCTION Icefront::CreatePVectorDiagnosticStokes{{{1*/
    625 void Icefront::CreatePVectorDiagnosticStokes( Vec pg){
     625ElementVector* Icefront::CreatePVectorDiagnosticStokes(void){
    626626
    627627        /*Constants*/
    … …  
    636636        double      xyz_list[NUMVERTICESQUA][3];
    637637        double      normal[3];
    638         double      pe_g[numdofs] = {0.0};
    639638        double      l1l4[4];
    640         int        *doflist = NULL;
    641639        GaussPenta *gauss = NULL;
    642         Penta      *penta = NULL;
    643 
    644         /* Get dof list and node coordinates: */
    645         penta=(Penta*)element;
    646         GetDofList(&doflist,StokesApproximationEnum,GsetEnum);
     640
     641        Penta* penta=(Penta*)element;
     642
     643        /*Initialize Element vector and return if necessary*/
     644        if(penta->IsOnWater()) return NULL;
     645        ElementVector* pe=NewElementVector(nodes,NUMVERTICESQUA,this->parameters);
     646
     647        /*Retrieve all inputs and parameters*/
    647648        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICESQUA);
    648        
    649         /*Recover inputs and parameters */
    650649        inputs->GetParameterValue(&fill,FillEnum);
    651650        rho_water=matpar->GetRhoWater();
    … …  
    687686                for(i=0;i<NUMVERTICESQUA;i++){
    688687                        for(j=0;j<NDOF4;j++){
    689                                 if(j<3)  pe_g[i*NDOF4+j]+=Jdet*gauss->weight*pressure*l1l4[i]*normal[j];
    690                                 else     pe_g[i*NDOF4+j]+=0; //pressure term
     688                                if(j<3)  pe->values[i*NDOF4+j]+=Jdet*gauss->weight*pressure*l1l4[i]*normal[j];
     689                                else     pe->values[i*NDOF4+j]+=0; //pressure term
    691690                        }
    692691                }
    693692        }
    694693
    695         VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES);
    696 
    697         /*Free ressources:*/
     694        /*Clean up and return*/
    698695        delete gauss;
    699         xfree((void**)&doflist);
     696        return pe;
    700697}
    701698/*}}}*/

  • issm/trunk/src/c/objects/Loads/Icefront.h

    r5909 r5935  
    7777                /*}}}*/
    7878                /*Load management: {{{1*/
    79                 void  CreatePVectorDiagnosticMacAyeal2d(Vec pg,Vec pf);
    80                 void  CreatePVectorDiagnosticMacAyeal3d(Vec pg);
    81                 void  CreatePVectorDiagnosticPattyn(Vec pg);
    82                 void  CreatePVectorDiagnosticStokes(Vec pg);
     79                ElementVector* CreatePVectorAdjointHoriz(void);
     80                ElementVector* CreatePVectorDiagnosticHoriz(void);
     81                ElementVector* CreatePVectorDiagnosticMacAyeal2d(void);
     82                ElementVector* CreatePVectorDiagnosticMacAyeal3d(void);
     83                ElementVector* CreatePVectorDiagnosticPattyn(void);
     84                ElementVector* CreatePVectorDiagnosticStokes(void);
    8385                void  GetDofList(int** pdoflist,int approximation_enum,int setenum);
    8486                void GetSegmentNormal(double* normal,double xyz_list[2][3]);
Note: See TracChangeset for help on using the changeset viewer.