Changeset 5926


Ignore:
Timestamp:
09/21/10 15:52:55 (14 years ago)
Author:
Mathieu Morlighem
Message:

more ElementVector

Location:
issm/trunk/src/c/objects
Files:
6 edited

Legend:

Unmodified
Added
Removed

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

    r5925 r5926  
    740740
    741741        /*retrive parameters: */
     742        ElementVector* pe=NULL;
    742743        int analysis_type;
    743744        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    … …  
    750751        switch(analysis_type){
    751752                case DiagnosticHorizAnalysisEnum:
    752                         CreatePVectorDiagnosticHoriz( pg);
     753                        pe=CreatePVectorDiagnosticHoriz();
     754                        if(pe) pe->AddToGlobal(pg,pf);
     755                        delete pe;
    753756                        break;
    754757                case AdjointHorizAnalysisEnum:
    755                         CreatePVectorAdjointHoriz( pg);
     758                        pe=CreatePVectorAdjointHoriz();
     759                        if(pe) pe->AddToGlobal(pg,pf);
     760                        delete pe;
    756761                        break;
    757762                case DiagnosticHutterAnalysisEnum:
    … …  
    29452950/*}}}*/
    29462951/*FUNCTION Penta::CreatePVectorAdjointHoriz{{{1*/
    2947 void  Penta::CreatePVectorAdjointHoriz(Vec pg){
     2952ElementVector* Penta::CreatePVectorAdjointHoriz(void){
    29482953
    29492954        int approximation;
    … …  
    29522957        switch(approximation){
    29532958                case MacAyealApproximationEnum:
    2954                         CreatePVectorAdjointMacAyeal( pg);
    2955                         break;
     2959                        return CreatePVectorAdjointMacAyeal();
    29562960                case PattynApproximationEnum:
    2957                         CreatePVectorAdjointPattyn( pg);
    2958                         break;
     2961                        return CreatePVectorAdjointPattyn();
    29592962                case NoneApproximationEnum:
    2960                         return;
     2963                        return NULL;
    29612964                case StokesApproximationEnum:
    2962                         CreatePVectorAdjointStokes( pg);
    2963                         break;
     2965                        return CreatePVectorAdjointStokes();
    29642966                default:
    29652967                        ISSMERROR("Approximation %s not supported yet",EnumToString(approximation));
    … …  
    29682970/*}}}*/
    29692971/*FUNCTION Penta::CreatePVectorAdjointMacAyeal{{{1*/
    2970 void  Penta::CreatePVectorAdjointMacAyeal(Vec p_g){
    2971 
    2972         if (!IsOnBed() || IsOnWater()) return;
     2972ElementVector* Penta::CreatePVectorAdjointMacAyeal(){
     2973
     2974        if (!IsOnBed() || IsOnWater()) return NULL;
    29732975
    29742976        /*Call Tria function*/
    … …  
    29762978        ElementVector* pe=tria->CreatePVectorAdjointHoriz();
    29772979        delete tria->matice; delete tria;
    2978         pe->AddToGlobal(p_g,NULL);
    2979         delete pe;
    29802980
    29812981        /*clean up and return*/
    2982         return;
     2982        return pe;
    29832983}
    29842984/*}}}*/
    29852985/*FUNCTION Penta::CreatePVectorAdjointPattyn{{{1*/
    2986 void  Penta::CreatePVectorAdjointPattyn(Vec p_g){
    2987 
    2988         if (!IsOnSurface() || IsOnWater()) return;
     2986ElementVector* Penta::CreatePVectorAdjointPattyn(void){
     2987
     2988        if (!IsOnSurface() || IsOnWater()) return NULL;
    29892989
    29902990        /*Call Tria function*/
    … …  
    29922992        ElementVector* pe=tria->CreatePVectorAdjointHoriz();
    29932993        delete tria->matice; delete tria;
    2994         pe->AddToGlobal(p_g,NULL);
    2995         delete pe;
    29962994
    29972995        /*clean up and return*/
    2998         return;
     2996        return pe;
    29992997}
    30002998/*}}}*/
    … …  
    30483046/*}}}*/
    30493047/*FUNCTION Penta::CreatePVectorCouplingPattynStokes {{{1*/
    3050 void Penta::CreatePVectorCouplingPattynStokes( Vec pg){
    3051 
    3052         /*indexing: */
     3048ElementVector* Penta::CreatePVectorCouplingPattynStokes(void){
     3049
     3050        /*compute all load vectorsfor this element*/
     3051        ElementVector* pe1=CreatePVectorCouplingPattynStokesViscous();
     3052        ElementVector* pe2=CreatePVectorCouplingPattynStokesFriction();
     3053        ElementVector* pe =new ElementVector(pe1,pe2);
     3054
     3055        /*clean-up and return*/
     3056        delete pe1;
     3057        delete pe2;
     3058        return pe;
     3059}
     3060/*}}}*/
     3061/*FUNCTION Penta::CreatePVectorCouplingPattynStokesViscous {{{1*/
     3062ElementVector* Penta::CreatePVectorCouplingPattynStokesViscous(void){
     3063
     3064        const int numdof=NUMVERTICES*NDOF4;
    30533065        int i,j;
    3054         const int numdof=NUMVERTICES*NDOF4;
    3055         int*      doflist=NULL;
    3056 
    3057         /*parameters: */
     3066        int     ig;
    30583067        double             xyz_list_tria[NUMVERTICES2D][3];
    30593068        double         xyz_list[NUMVERTICES][3];
    30603069        double             bed_normal[3];
    3061 
    3062         /* gaussian points: */
    3063         int     ig;
    30643070        GaussPenta *gauss=NULL;
    3065 
    30663071        double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    30673072        double  viscosity, w, alpha2_gauss;
    30683073        double  dw[3];
    3069 
    3070         /*matrices: */
    30713074        double     Pe_gaussian[24]={0.0}; //for the six nodes
    30723075        double     l1l6[6]; //for the six nodes of the penta
    … …  
    30743077        double     Jdet;
    30753078        double     Jdet2d;
    3076 
    30773079        Tria*     tria=NULL;
    30783080        Friction* friction=NULL;
    3079 
    3080         /*parameters: */
    30813081        double stokesreconditioning;
    30823082        int    approximation;
    30833083        int    analysis_type;
    30843084
    3085         /*retrieve inputs :*/
     3085        /*Initialize Element vector and return if necessary*/
     3086        if(IsOnWater()) return NULL;
    30863087        inputs->GetParameterValue(&approximation,ApproximationEnum);
    3087 
    3088         /*retrieve some parameters: */
     3088        if(approximation!=PattynStokesApproximationEnum) return NULL;
     3089        ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     3090
     3091        /*Retrieve all inputs and parameters*/
     3092        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    30893093        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    30903094        this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
    3091 
    3092         /*If on water or not Stokes, skip load: */
    3093         if(IsOnWater() || approximation!=PattynStokesApproximationEnum) return;
    3094 
    3095         /* Get node coordinates and dof list: */
    3096         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3097         GetDofList(&doflist,StokesApproximationEnum,GsetEnum);
    3098 
    3099         /*Retrieve all inputs we will be needing: */
    31003095        Input* vx_input=inputs->GetInput(VxEnum);               ISSMASSERT(vx_input);
    31013096        Input* vy_input=inputs->GetInput(VyEnum);               ISSMASSERT(vy_input);
    … …  
    31093104                gauss->GaussPoint(ig);
    31103105
    3111                 /*Compute strain rate and viscosity: */
    31123106                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    31133107                matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
    31143108
    3115                 /* Get Jacobian determinant: */
    31163109                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    31173110
    3118                 /* Get nodal functions */
    31193111                GetNodalFunctionsP1(&l1l6[0], gauss);
    31203112                GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
    31213113
    3122                 /*Compute the derivative of VzPattyn*/
    31233114                vzpattyn_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
    31243115
    3125                 /* Build gaussian vector */
    31263116                for(i=0;i<NUMVERTICES;i++){
    3127                         Pe_gaussian[i*NDOF4+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dh1dh6[2][i]/2;
    3128                         Pe_gaussian[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dh1dh6[2][i]/2;
    3129                         Pe_gaussian[i*NDOF4+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dh1dh6[0][i]+dw[1]*dh1dh6[1][i]+dw[2]*dh1dh6[2][i])/2;
    3130                         Pe_gaussian[i*NDOF4+3]+=Jdet*gauss->weight*stokesreconditioning*dw[2]*l1l6[i];
    3131                 }
    3132         }
     3117                        pe->values[i*NDOF4+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dh1dh6[2][i]/2;
     3118                        pe->values[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dh1dh6[2][i]/2;
     3119                        pe->values[i*NDOF4+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dh1dh6[0][i]+dw[1]*dh1dh6[1][i]+dw[2]*dh1dh6[2][i])/2;
     3120                        pe->values[i*NDOF4+3]+=Jdet*gauss->weight*stokesreconditioning*dw[2]*l1l6[i];
     3121                }
     3122        }
     3123
     3124        /*Clean up and return*/
    31333125        delete gauss;
    3134 
    3135         /*Deal with 2d friction at the bedrock interface: */
    3136         if(IsOnBed() && !IsOnShelf()){
     3126        return pe;
     3127}
     3128/*}}}*/
     3129/*FUNCTION Penta::CreatePVectorCouplingPattynStokesFriction{{{1*/
     3130ElementVector* Penta::CreatePVectorCouplingPattynStokesFriction(void){
     3131
     3132        const int numdof=NUMVERTICES*NDOF4;
     3133        int i,j;
     3134        int     ig;
     3135        double             xyz_list_tria[NUMVERTICES2D][3];
     3136        double         xyz_list[NUMVERTICES][3];
     3137        double             bed_normal[3];
     3138        GaussPenta *gauss=NULL;
     3139        double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
     3140        double  viscosity, w, alpha2_gauss;
     3141        double  dw[3];
     3142        double     Pe_gaussian[24]={0.0}; //for the six nodes
     3143        double     l1l6[6]; //for the six nodes of the penta
     3144        double     dh1dh6[3][6]; //for the six nodes of the penta
     3145        double     Jdet;
     3146        double     Jdet2d;
     3147        Tria*     tria=NULL;
     3148        Friction* friction=NULL;
     3149        double stokesreconditioning;
     3150        int    approximation;
     3151        int    analysis_type;
     3152
     3153        /*Initialize Element vector and return if necessary*/
     3154        if(IsOnWater() || !IsOnBed() || IsOnShelf()) return NULL;
     3155        inputs->GetParameterValue(&approximation,ApproximationEnum);
     3156        if(approximation!=PattynStokesApproximationEnum) return NULL;
     3157        ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     3158
     3159        /*Retrieve all inputs and parameters*/
     3160        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     3161        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     3162        this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
     3163        Input* vx_input=inputs->GetInput(VxEnum);               ISSMASSERT(vx_input);
     3164        Input* vy_input=inputs->GetInput(VyEnum);               ISSMASSERT(vy_input);
     3165        Input* vz_input=inputs->GetInput(VzEnum);               ISSMASSERT(vz_input);
     3166        Input* vzpattyn_input=inputs->GetInput(VzPattynEnum);   ISSMASSERT(vzpattyn_input);
     3167
     3168
     3169        for(i=0;i<NUMVERTICES2D;i++){
     3170                for(j=0;j<3;j++){
     3171                        xyz_list_tria[i][j]=xyz_list[i][j];
     3172                }
     3173        }
     3174
     3175        /*build friction object, used later on: */
     3176        friction=new Friction("3d",inputs,matpar,analysis_type);
     3177
     3178        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
     3179        gauss=new GaussPenta(0,1,2,2);
     3180        for(ig=gauss->begin();ig<gauss->end();ig++){
     3181
     3182                gauss->GaussPoint(ig);
     3183
     3184                /*Get the Jacobian determinant */
     3185                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
     3186
     3187                /*Get L matrix : */
     3188                GetNodalFunctionsP1(l1l6, gauss);
     3189
     3190                /*Get normal vecyor to the bed */
     3191                BedNormal(&bed_normal[0],xyz_list_tria);
     3192
     3193                /*Get Viscosity*/
     3194                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     3195                matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     3196
     3197                /*Get friction*/
     3198                friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum);
     3199
     3200                /*Get w and its derivatives*/
     3201                vzpattyn_input->GetParameterValue(&w, gauss);
     3202                vzpattyn_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
    31373203
    31383204                for(i=0;i<NUMVERTICES2D;i++){
    3139                         for(j=0;j<3;j++){
    3140                                 xyz_list_tria[i][j]=xyz_list[i][j];
    3141                         }
    3142                 }
    3143 
    3144                 /*build friction object, used later on: */
    3145                 friction=new Friction("3d",inputs,matpar,analysis_type);
    3146 
    3147                 /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    3148                 gauss=new GaussPenta(0,1,2,2);
    3149                 for(ig=gauss->begin();ig<gauss->end();ig++){
    3150 
    3151                         gauss->GaussPoint(ig);
    3152 
    3153                         /*Get the Jacobian determinant */
    3154                         GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    3155 
    3156                         /*Get L matrix : */
    3157                         GetNodalFunctionsP1(l1l6, gauss);
    3158 
    3159                         /*Get normal vecyor to the bed */
    3160                         BedNormal(&bed_normal[0],xyz_list_tria);
    3161 
    3162                         /*Get Viscosity*/
    3163                         this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    3164                         matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
    3165 
    3166                         /*Get friction*/
    3167                         friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum);
    3168 
    3169                         /*Get w and its derivatives*/
    3170                         vzpattyn_input->GetParameterValue(&w, gauss);
    3171                         vzpattyn_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
    3172 
    3173                         for(i=0;i<NUMVERTICES2D;i++){
    3174                                 Pe_gaussian[i*NDOF4+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+viscosity*dw[2]*bed_normal[0])*l1l6[i];
    3175                                 Pe_gaussian[i*NDOF4+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+viscosity*dw[2]*bed_normal[1])*l1l6[i];
    3176                                 Pe_gaussian[i*NDOF4+2]+=Jdet2d*gauss->weight*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*l1l6[i];
    3177                         }
    3178                 }
    3179                 /*Free ressources:*/
    3180                 delete gauss;
    3181         } //if ( (IsOnBed()==1) && (IsOnShelf()==1))
    3182 
    3183         /*Add Pe_reduced to global vector pg: */
    3184         VecSetValues(pg,numdof,doflist,(const double*)Pe_gaussian,ADD_VALUES);
    3185 
    3186         /*Free ressources:*/
    3187         xfree((void**)&doflist);
    3188 
     3205                        pe->values[i*NDOF4+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+viscosity*dw[2]*bed_normal[0])*l1l6[i];
     3206                        pe->values[i*NDOF4+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+viscosity*dw[2]*bed_normal[1])*l1l6[i];
     3207                        pe->values[i*NDOF4+2]+=Jdet2d*gauss->weight*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*l1l6[i];
     3208                }
     3209        }
     3210
     3211        /*Clean up and return*/
     3212        delete gauss;
     3213        return pe;
    31893214}
    31903215/*}}}*/
    31913216/*FUNCTION Penta::CreatePVectorDiagnosticHoriz{{{1*/
    3192 void  Penta::CreatePVectorDiagnosticHoriz(Vec pg){
     3217ElementVector* Penta::CreatePVectorDiagnosticHoriz(void){
    31933218
    31943219        int approximation;
    … …  
    31973222        switch(approximation){
    31983223                case MacAyealApproximationEnum:
    3199                         CreatePVectorDiagnosticMacAyeal( pg);
    3200                         break;
     3224                        return CreatePVectorDiagnosticMacAyeal();
    32013225                case PattynApproximationEnum:
    3202                         CreatePVectorDiagnosticPattyn( pg);
    3203                         break;
     3226                        return CreatePVectorDiagnosticPattyn();
    32043227                case HutterApproximationEnum:
    3205                         return;
     3228                        return NULL;
    32063229                case NoneApproximationEnum:
    3207                         return;
     3230                        return NULL;
    32083231                case StokesApproximationEnum:
    3209                         CreatePVectorDiagnosticStokes( pg);
    3210                         break;
     3232                        return CreatePVectorDiagnosticStokes();
    32113233                case MacAyealPattynApproximationEnum:
    3212                         CreatePVectorDiagnosticMacAyeal( pg);
    3213                         CreatePVectorDiagnosticPattyn( pg);
    3214                         break;
     3234                        return CreatePVectorDiagnosticMacAyealPattyn();
    32153235                case PattynStokesApproximationEnum:
    3216                         CreatePVectorDiagnosticPattyn( pg);
    3217                         CreatePVectorDiagnosticStokes( pg);
    3218                         CreatePVectorCouplingPattynStokes( pg);
    3219                         break;
     3236                        return CreatePVectorDiagnosticPattynStokes();
    32203237                default:
    32213238                        ISSMERROR("Approximation %s not supported yet",EnumToString(approximation));
    32223239        }
     3240}
     3241/*}}}*/
     3242/*FUNCTION Penta::CreatePVectorDiagnosticMacAyealPattyn{{{1*/
     3243ElementVector* Penta::CreatePVectorDiagnosticMacAyealPattyn(void){
     3244
     3245        /*compute all load vectorsfor this element*/
     3246        ElementVector* pe1=CreatePVectorDiagnosticMacAyeal();
     3247        ElementVector* pe2=CreatePVectorDiagnosticPattyn();
     3248        ElementVector* pe =new ElementVector(pe1,pe2);
     3249
     3250        /*clean-up and return*/
     3251        delete pe1;
     3252        delete pe2;
     3253        return pe;
     3254}
     3255/*}}}*/
     3256/*FUNCTION Penta::CreatePVectorDiagnosticPattynStokes{{{1*/
     3257ElementVector* Penta::CreatePVectorDiagnosticPattynStokes(void){
     3258
     3259        /*compute all load vectorsfor this element*/
     3260        ElementVector* pe1=CreatePVectorDiagnosticPattyn();
     3261        ElementVector* pe2=CreatePVectorDiagnosticStokes();
     3262        ElementVector* pe3=CreatePVectorCouplingPattynStokes();
     3263        ElementVector* pe =new ElementVector(pe1,pe2,pe3);
     3264
     3265        /*clean-up and return*/
     3266        delete pe1;
     3267        delete pe2;
     3268        delete pe3;
     3269        return pe;
    32233270}
    32243271/*}}}*/
    … …  
    33293376/*}}}*/
    33303377/*FUNCTION Penta::CreatePVectorDiagnosticMacAyeal{{{1*/
    3331 void  Penta::CreatePVectorDiagnosticMacAyeal(Vec pg){
    3332 
    3333         /*Figure out if this penta is collapsed. If so, then bailout, except if it is at the
    3334           bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build
    3335           the stiffness matrix. */
    3336         if (!IsOnBed() || IsOnWater()) return;
     3378ElementVector* Penta::CreatePVectorDiagnosticMacAyeal(void){
     3379
     3380        if (!IsOnBed() || IsOnWater()) return NULL;
    33373381
    33383382        /*Call Tria function*/
    … …  
    33403384        ElementVector* pe=tria->CreatePVectorDiagnosticMacAyeal();
    33413385        delete tria->matice; delete tria;
    3342         pe->AddToGlobal(pg,NULL);
    3343         delete pe;
    3344 
    3345         /*clean up and return*/
    3346         return;
     3386
     3387        /*Clean up and return*/
     3388        return pe;
    33473389}
    33483390/*}}}*/
    33493391/*FUNCTION Penta::CreatePVectorDiagnosticPattyn{{{1*/
    3350 void Penta::CreatePVectorDiagnosticPattyn( Vec pg){
    3351 
    3352         int i,j;
     3392ElementVector* Penta::CreatePVectorDiagnosticPattyn(void){
    33533393
    33543394        /* node data: */
    33553395        const int    numdof=NDOF2*NUMVERTICES;
     3396        int i,j;
     3397        int     ig;
    33563398        double       xyz_list[NUMVERTICES][3];
    3357         int*         doflist=NULL;
    3358 
    3359         /* parameters: */
    33603399        double  slope[3]; //do not put 2! this goes into GetParameterDerivativeValue, which addresses slope[3] also!
    33613400        double  driving_stress_baseline;
    33623401        double  thickness;
    3363 
    3364         /* gaussian points: */
    3365         int     ig;
    33663402        GaussPenta *gauss=NULL;
    3367 
    3368         /* Jacobian: */
    33693403        double Jdet;
    3370 
    3371         /*nodal functions: */
    33723404        double l1l6[6];
    3373 
    3374         /*element vector at the gaussian points: */
    3375         double  pe_g[numdof]={0.0};
    33763405        double  pe_g_gaussian[numdof];
    33773406
    3378         /*If on water, skip load: */
    3379         if(IsOnWater())return;
    3380 
    3381         /*Implement standard penta element: */
    3382 
    3383         /* Get node coordinates and dof list: */
     3407        /*Initialize Element vector and return if necessary*/
     3408        if(IsOnWater()) return NULL;
     3409        ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters,PattynApproximationEnum);
     3410
     3411        /*Retrieve all inputs and parameters*/
    33843412        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3385         GetDofList(&doflist,PattynApproximationEnum,GsetEnum);
    3386 
    3387         /*Retrieve all inputs we will be needing: */
    33883413        Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input);
    33893414        Input* surface_input=inputs->GetInput(SurfaceEnum);     ISSMASSERT(surface_input);
    … …  
    33953420                gauss->GaussPoint(ig);
    33963421
    3397                 /*Compute thickness at gaussian point: */
    33983422                thickness_input->GetParameterValue(&thickness, gauss);
    3399 
    3400                 /*Compute slope at gaussian point: */
    34013423                surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
    34023424
    3403                 /* Get Jacobian determinant: */
    34043425                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3405 
    3406                 /*Get nodal functions: */
    34073426                GetNodalFunctionsP1(l1l6, gauss);
    34083427
    3409                 /*Compute driving stress: */
    34103428                driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG();
    34113429
    3412                 /*Build pe_g_gaussian vector: */
    3413                 for (i=0;i<NUMVERTICES;i++){
    3414                         for (j=0;j<NDOF2;j++){
    3415                                 pe_g_gaussian[i*NDOF2+j]=-driving_stress_baseline*slope[j]*Jdet*gauss->weight*l1l6[i];
    3416                         }
    3417                 }
    3418 
    3419                 /*Add pe_g_gaussian vector to pe_g: */
    3420                 for( i=0; i<numdof; i++)pe_g[i]+=pe_g_gaussian[i];
    3421         }
    3422 
    3423         /*Add pe_g to global vector pg: */
    3424         VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    3425 
    3426         xfree((void**)&doflist);
     3430                for(i=0;i<NUMVERTICES;i++) for(j=0;j<NDOF2;j++) pe->values[i*NDOF2+j]+= -driving_stress_baseline*slope[j]*Jdet*gauss->weight*l1l6[i];
     3431        }
     3432
     3433        /*Clean up and return*/
    34273434        delete gauss;
     3435        return pe;
    34283436}
    34293437/*}}}*/
    34303438/*FUNCTION Penta::CreatePVectorDiagnosticStokes {{{1*/
    3431 void Penta::CreatePVectorDiagnosticStokes( Vec pg){
    3432 
    3433         /*indexing: */
     3439ElementVector* Penta::CreatePVectorDiagnosticStokes(void){
     3440
     3441        /*compute all load vectorsfor this element*/
     3442        ElementVector* pe1=CreatePVectorDiagnosticStokesViscous();
     3443        ElementVector* pe2=CreatePVectorDiagnosticStokesShelf();
     3444        ElementVector* pe =new ElementVector(pe1,pe2);
     3445
     3446        /*clean-up and return*/
     3447        delete pe1;
     3448        delete pe2;
     3449        return pe;
     3450}
     3451/*}}}*/
     3452/*FUNCTION Penta::CreatePVectorDiagnosticStokesViscous {{{1*/
     3453ElementVector* Penta::CreatePVectorDiagnosticStokesViscous(void){
     3454
    34343455        int i,j;
     3456        int     ig;
    34353457        const int numdof=NUMVERTICES*NDOF4;
    34363458        int       numdof2d=NUMVERTICES2D*NDOF4;
    3437         int*      doflist=NULL;
    3438 
    3439         /*Material properties: */
    34403459        double         gravity,rho_ice,rho_water;
    3441 
    3442         /*parameters: */
    34433460        double             xyz_list_tria[NUMVERTICES2D][3];
    34443461        double         xyz_list[NUMVERTICES][3];
    34453462        double             bed_normal[3];
    34463463        double         bed;
    3447 
    3448         /* gaussian points: */
    3449         int     ig;
    34503464        GaussPenta *gauss=NULL;
    3451 
    34523465        double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    34533466        double  viscosity;
    34543467        double  water_pressure;
    3455 
    3456         /*matrices: */
    34573468        double     Pe_temp[27]={0.0}; //for the six nodes and the bubble
    34583469        double     Pe_gaussian[27]={0.0}; //for the six nodes and the bubble
    … …  
    34703481        double     D_scalar;
    34713482        double     tBD[27][8];
    3472 
    34733483        Tria*            tria=NULL;
    3474 
    3475         /*parameters: */
    34763484        double stokesreconditioning;
    3477 
    3478         /*inputs: */
    34793485        int  approximation;
    34803486
    3481         /*retrieve inputs :*/
     3487        /*Initialize Element vector and return if necessary*/
     3488        if(IsOnWater()) return NULL;
    34823489        inputs->GetParameterValue(&approximation,ApproximationEnum);
    3483 
    3484         /*retrieve some parameters: */
     3490        if(approximation!=StokesApproximationEnum) return NULL;
     3491        ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     3492
     3493        /*Retrieve all inputs and parameters*/
    34853494        this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
    3486 
    3487         /*If on water or not Stokes, skip load: */
    3488         if(IsOnWater() || approximation!=StokesApproximationEnum) return;
    3489 
    3490         /*recovre material parameters: */
    34913495        rho_water=matpar->GetRhoWater();
    34923496        rho_ice=matpar->GetRhoIce();
    34933497        gravity=matpar->GetG();
    3494 
    3495         /* Get node coordinates and dof list: */
    34963498        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3497         GetDofList(&doflist,StokesApproximationEnum,GsetEnum);
    3498 
    3499         /*Retrieve all inputs we will be needing: */
    35003499        Input* vx_input=inputs->GetInput(VxEnum);   ISSMASSERT(vx_input);
    35013500        Input* vy_input=inputs->GetInput(VyEnum);   ISSMASSERT(vy_input);
    … …  
    35093508                gauss->GaussPoint(ig);
    35103509
    3511                 /*Compute strain rate and viscosity: */
    35123510                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    35133511                matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
    3514 
    3515                 /* Get Jacobian determinant: */
    35163512                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3517 
    3518                 /* Get nodal functions */
    35193513                GetNodalFunctionsMINI(&l1l7[0], gauss);
    35203514
    … …  
    35473541                for(i=0;i<27;i++) for(j=0;j<3;j++) Ke_temp[i][j]+=Ke_gaussian[i][j];
    35483542        }
     3543
     3544        /*Condensation*/
     3545        ReduceVectorStokes(pe->values, &Ke_temp[0][0], &Pe_temp[0]);
     3546
     3547        /*Clean up and return*/
    35493548        delete gauss;
    3550 
    3551         /*Deal with 2d friction at the bedrock interface: */
    3552         if(IsOnBed() && IsOnShelf()){
    3553 
    3554                 for(i=0;i<NUMVERTICES2D;i++){
    3555                         for(j=0;j<3;j++){
    3556                                 xyz_list_tria[i][j]=xyz_list[i][j];
    3557                         }
    3558                 }
    3559 
    3560                 /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    3561                 gauss=new GaussPenta(0,1,2,2);
    3562                 for(ig=gauss->begin();ig<gauss->end();ig++){
    3563 
    3564                         gauss->GaussPoint(ig);
    3565 
    3566                         /*Get the Jacobian determinant */
    3567                         GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    3568 
    3569                         /* Get bed at gaussian point */
    3570                         bed_input->GetParameterValue(&bed, gauss);
    3571 
    3572                         /*Get L matrix : */
    3573                         GetNodalFunctionsP1(l1l6, gauss);
    3574 
    3575                         /*Get water_pressure at gaussian point */
    3576                         water_pressure=gravity*rho_water*bed;
    3577 
    3578                         /*Get normal vecyor to the bed */
    3579                         BedNormal(&bed_normal[0],xyz_list_tria);
    3580 
    3581                         for(i=0;i<NUMVERTICES2D;i++){
    3582                                 for(j=0;j<3;j++){
    3583                                         Pe_temp[i*NDOF4+j]+=water_pressure*gauss->weight*Jdet2d*l1l6[i]*bed_normal[j];
    3584                                 }
    3585                         }
    3586                 }
    3587                 /*Free ressources:*/
    3588                 delete gauss;
    3589         } //if ( (IsOnBed()==1) && (IsOnShelf()==1))
    3590 
    3591         /*Reduce the matrix */
    3592         ReduceVectorStokes(&Pe_reduced[0], &Ke_temp[0][0], &Pe_temp[0]);
    3593 
    3594         /*Add Pe_reduced to global vector pg: */
    3595         VecSetValues(pg,numdof,doflist,(const double*)Pe_reduced,ADD_VALUES);
    3596 
    3597         /*Free ressources:*/
    3598         xfree((void**)&doflist);
    3599 
     3549        return pe;
     3550}
     3551/*}}}*/
     3552/*FUNCTION Penta::CreatePVectorDiagnosticStokesShelf{{{1*/
     3553ElementVector* Penta::CreatePVectorDiagnosticStokesShelf(void){
     3554
     3555        int i,j;
     3556        int     ig;
     3557        const int numdof=NUMVERTICES*NDOF4;
     3558        int       numdof2d=NUMVERTICES2D*NDOF4;
     3559        double         gravity,rho_ice,rho_water;
     3560        double             xyz_list_tria[NUMVERTICES2D][3];
     3561        double         xyz_list[NUMVERTICES][3];
     3562        double             bed_normal[3];
     3563        double         bed;
     3564        GaussPenta *gauss=NULL;
     3565        double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
     3566        double  viscosity;
     3567        double  water_pressure;
     3568        double     Pe_temp[27]={0.0}; //for the six nodes and the bubble
     3569        double     Pe_gaussian[27]={0.0}; //for the six nodes and the bubble
     3570        double     Ke_temp[27][3]={0.0}; //for the six nodes and the bubble
     3571        double     Pe_reduced[numdof]; //for the six nodes only
     3572        double     Ke_gaussian[27][3];
     3573        double     l1l6[6]; //for the six nodes of the penta
     3574        double     l1l7[7]; //for the six nodes and the bubble
     3575        double     B[8][27];
     3576        double     B_prime[8][27];
     3577        double     B_prime_bubble[8][3];
     3578        double     Jdet;
     3579        double     Jdet2d;
     3580        double     D[8][8]={0.0};
     3581        double     D_scalar;
     3582        double     tBD[27][8];
     3583        Tria*            tria=NULL;
     3584        double stokesreconditioning;
     3585        int  approximation;
     3586
     3587        /*Initialize Element vector and return if necessary*/
     3588        if(IsOnWater() || !IsOnBed() || !IsOnShelf()) return NULL;
     3589        inputs->GetParameterValue(&approximation,ApproximationEnum);
     3590        if(approximation!=StokesApproximationEnum) return NULL;
     3591        ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     3592
     3593        /*Retrieve all inputs and parameters*/
     3594        this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
     3595        rho_water=matpar->GetRhoWater();
     3596        rho_ice=matpar->GetRhoIce();
     3597        gravity=matpar->GetG();
     3598        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     3599        Input* vx_input=inputs->GetInput(VxEnum);   ISSMASSERT(vx_input);
     3600        Input* vy_input=inputs->GetInput(VyEnum);   ISSMASSERT(vy_input);
     3601        Input* vz_input=inputs->GetInput(VzEnum);   ISSMASSERT(vz_input);
     3602        Input* bed_input=inputs->GetInput(BedEnum); ISSMASSERT(bed_input);
     3603
     3604
     3605        for(i=0;i<NUMVERTICES2D;i++){
     3606                for(j=0;j<3;j++){
     3607                        xyz_list_tria[i][j]=xyz_list[i][j];
     3608                }
     3609        }
     3610
     3611        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
     3612        gauss=new GaussPenta(0,1,2,2);
     3613        for(ig=gauss->begin();ig<gauss->end();ig++){
     3614
     3615                gauss->GaussPoint(ig);
     3616
     3617                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
     3618                bed_input->GetParameterValue(&bed, gauss);
     3619                GetNodalFunctionsP1(l1l6, gauss);
     3620
     3621                water_pressure=gravity*rho_water*bed;
     3622
     3623                BedNormal(&bed_normal[0],xyz_list_tria);
     3624
     3625                for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) pe->values[i*NDOF4+j]+=water_pressure*gauss->weight*Jdet2d*l1l6[i]*bed_normal[j];
     3626        }
     3627
     3628        /*Clean up and return*/
     3629        delete gauss;
     3630        return pe;
    36003631}
    36013632/*}}}*/
    36023633/*FUNCTION Penta::CreatePVectorAdjointStokes{{{1*/
    3603 void  Penta::CreatePVectorAdjointStokes(Vec p_g){
    3604 
    3605         Tria* tria=NULL;
    3606 
    3607         /*If on water, skip: */
    3608         if(IsOnWater() || !IsOnSurface())return;
    3609 
    3610         /*Call Tria's function*/
    3611         tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
    3612         tria->CreatePVectorAdjointStokes(p_g);
     3634ElementVector* Penta::CreatePVectorAdjointStokes(void){
     3635
     3636        if (!IsOnSurface() || IsOnWater()) return NULL;
     3637
     3638        /*Call Tria function*/
     3639        Tria* tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
     3640        ElementVector* pe=tria->CreatePVectorAdjointStokes();
    36133641        delete tria->matice; delete tria;
    3614         return;
     3642
     3643        /*clean up and return*/
     3644        return pe;
    36153645}
    36163646/*}}}*/

  • TabularUnified issm/trunk/src/c/objects/Elements/Penta.h

    r5917 r5926  
    152152                void      CreatePVectorBalancedthickness( Vec pg);
    153153                void      CreatePVectorBalancedvelocities( Vec pg);
    154                 void      CreatePVectorAdjointHoriz( Vec pg);
    155                 void      CreatePVectorAdjointMacAyeal( Vec pg);
    156                 void      CreatePVectorAdjointPattyn( Vec pg);
    157                 void      CreatePVectorCouplingPattynStokes( Vec pg);
    158                 void      CreatePVectorDiagnosticHoriz( Vec pg);
     154                ElementVector* CreatePVectorAdjointHoriz(void);
     155                ElementVector* CreatePVectorAdjointMacAyeal(void);
     156                ElementVector* CreatePVectorAdjointPattyn(void);
     157                ElementVector* CreatePVectorCouplingPattynStokes(void);
     158                ElementVector* CreatePVectorCouplingPattynStokesViscous(void);
     159                ElementVector* CreatePVectorCouplingPattynStokesFriction(void);
     160                ElementVector* CreatePVectorDiagnosticHoriz(void);
    159161                void      CreatePVectorDiagnosticHutter( Vec pg);
    160                 void      CreatePVectorDiagnosticMacAyeal( Vec pg);
    161                 void      CreatePVectorDiagnosticPattyn( Vec pg);
    162                 void      CreatePVectorDiagnosticStokes( Vec pg);
    163                 void      CreatePVectorAdjointStokes( Vec pg);
     162                ElementVector* CreatePVectorDiagnosticMacAyeal(void);
     163                ElementVector* CreatePVectorDiagnosticMacAyealPattyn(void);
     164                ElementVector* CreatePVectorDiagnosticPattyn(void);
     165                ElementVector* CreatePVectorDiagnosticPattynStokes(void);
     166                ElementVector* CreatePVectorDiagnosticStokes(void);
     167                ElementVector* CreatePVectorDiagnosticStokesViscous(void);
     168                ElementVector* CreatePVectorDiagnosticStokesShelf(void);
     169                ElementVector* CreatePVectorAdjointStokes(void);
    164170                void      CreatePVectorDiagnosticVert( Vec pg);
    165171                void      CreatePVectorMelting( Vec pg);

  • TabularUnified issm/trunk/src/c/objects/Elements/Tria.cpp

    r5924 r5926  
    39683968/*}}}*/
    39693969/*FUNCTION Tria::CreatePVectorAdjointStokes{{{1*/
    3970 void Tria::CreatePVectorAdjointStokes(Vec p_g){
    3971 
     3970ElementVector* Tria::CreatePVectorAdjointStokes(void){
     3971
     3972        const int    numdof=NDOF4*NUMVERTICES;
    39723973        int i;
    3973 
    3974         /* node data: */
    3975         const int    numdof=NDOF4*NUMVERTICES;
     3974        int     ig;
    39763975        double       xyz_list[NUMVERTICES][3];
    3977         int*         doflist=NULL;
    3978 
    3979         /* grid data: */
    39803976        double vx_list[NUMVERTICES];
    39813977        double vy_list[NUMVERTICES];
    … …  
    39853981        double duy_list[NUMVERTICES];
    39863982        double weights_list[NUMVERTICES];
    3987 
    3988         /* gaussian points: */
    3989         int     ig;
    39903983        GaussTria* gauss=NULL;
    3991 
    3992         /* parameters: */
    39933984        double  obs_velocity_mag,velocity_mag;
    39943985        double  dux,duy;
    39953986        double  meanvel, epsvel;
    3996 
    3997         /*element vector : */
    39983987        double  pe_g[numdof]={0.0};
    3999 
    4000         /* Jacobian: */
    40013988        double Jdet;
    4002 
    4003         /*nodal functions: */
    40043989        double l1l2l3[3];
    4005 
    4006         /*relative and algorithmic fitting: */
    40073990        double scalex=0;
    40083991        double scaley=0;
    … …  
    40123995        int    response;
    40133996
    4014         /*retrieve some parameters: */
     3997        /*Initialize Element vector and return if necessary*/
     3998        if(IsOnWater()) return NULL;
     3999        ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     4000
     4001        /*Retrieve all inputs and parameters*/
     4002        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    40154003        this->parameters->FindParam(&meanvel,MeanVelEnum);
    40164004        this->parameters->FindParam(&epsvel,EpsVelEnum);
    4017 
    4018         /* Get node coordinates and dof list: */
    4019         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    4020         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    4021 
    4022         /* Recover input data: */
    40234005        GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum);
    40244006        GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum);
    … …  
    40264008        GetParameterListOnVertices(&vy_list[0],VyEnum);
    40274009        GetParameterListOnVertices(&weights_list[0],WeightsEnum);
    4028 
    40294010        inputs->GetParameterValue(&response,CmResponseEnum);
    40304011        if(response==SurfaceAverageVelMisfitEnum){
    … …  
    41564137                gauss->GaussPoint(ig);
    41574138
    4158                 /* Get Jacobian determinant: */
    41594139                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    4160 
    4161                 /* Get nodal functions value at gaussian point:*/
    41624140                GetNodalFunctions(l1l2l3, gauss);
    41634141
    4164                 /*Build due_g_gaussian vector: we have three cases here, according to which type of misfit we are using. */
    4165 
    4166                 /*Compute absolute(x/y) at gaussian point: */
    41674142                TriaRef::GetParameterValue(&dux, &dux_list[0],gauss);
    41684143                TriaRef::GetParameterValue(&duy, &duy_list[0],gauss);
    41694144
    4170                 /*compute Du*/
    41714145                for (i=0;i<NUMVERTICES;i++){
    4172                         pe_g[i*NDOF4+0]+=dux*Jdet*gauss->weight*l1l2l3[i];
    4173                         pe_g[i*NDOF4+1]+=duy*Jdet*gauss->weight*l1l2l3[i];
     4146                        pe->values[i*NDOF4+0]+=dux*Jdet*gauss->weight*l1l2l3[i];
     4147                        pe->values[i*NDOF4+1]+=duy*Jdet*gauss->weight*l1l2l3[i];
    41744148                }
    41754149        }
    4176 
    4177         /*Add pe_g to global vector p_g: */
    4178         VecSetValues(p_g,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    41794150
    41804151        /*Clean up and return*/
    41814152        delete gauss;
    4182         xfree((void**)&doflist);
     4153        return pe;
    41834154}
    41844155/*}}}*/

  • TabularUnified issm/trunk/src/c/objects/Elements/Tria.h

    r5924 r5926  
    145145                ElementVector* CreatePVectorDiagnosticMacAyeal(void);
    146146                ElementVector* CreatePVectorAdjointHoriz(void);
    147                 void      CreatePVectorAdjointStokes(Vec pg);
     147                ElementVector* CreatePVectorAdjointStokes(void);
    148148                ElementVector* CreatePVectorAdjointBalancedthickness(void);
    149149                ElementVector* CreatePVectorDiagnosticHutter(void);

  • TabularUnified issm/trunk/src/c/objects/Numerics/ElementVector.cpp

    r5827 r5926  
    3131}
    3232/*}}}*/
     33/*FUNCTION ElementVector::ElementVector(ElementVector* pe1, ElementVector* pe2){{{1*/
     34ElementVector::ElementVector(ElementVector* pe1, ElementVector* pe2){
     35
     36        /*intermediaries*/
     37        int i,j,counter;
     38        int gsize,fsize,ssize;
     39        int* P=NULL;
     40        bool found;
     41
     42        /*If one of the two matrix is NULL, we copy the other one*/
     43        if(!pe1 && !pe2){
     44                ISSMERROR("Two input element matrices are NULL");
     45        }
     46        else if(!pe1){
     47                this->Init(pe2);
     48                return;
     49        }
     50        else if(!pe2){
     51                this->Init(pe1);
     52                return;
     53        }
     54
     55        /*Initialize itransformation matrix pe[P[i]] = pe2[i]*/
     56        P=(int*)xmalloc(pe2->nrows*sizeof(int));
     57
     58        /*1: Get the new numbering of pe2 and get size of the new matrix*/
     59        gsize=pe1->nrows;
     60        for(i=0;i<pe2->nrows;i++){
     61                found=false;
     62                for(j=0;j<pe1->nrows;j++){
     63                        if(pe2->gglobaldoflist[i]==pe1->gglobaldoflist[j]){
     64                                found=true; P[i]=j; break;
     65                        }
     66                }
     67                if(!found){
     68                        P[i]=gsize; gsize++;
     69                }
     70        }
     71
     72        /*2: Initialize static fields*/
     73        this->nrows=gsize;
     74        this->pf=pe1->pf;
     75
     76        /*Gset and values*/
     77        this->gglobaldoflist=(int*)xmalloc(this->nrows*sizeof(int));
     78        this->values=(double*)xcalloc(this->nrows,sizeof(double));
     79        for(i=0;i<pe1->nrows;i++){
     80                this->values[i] += pe1->values[i];
     81                this->gglobaldoflist[i]=pe1->gglobaldoflist[i];
     82        }
     83        for(i=0;i<pe2->nrows;i++){
     84                this->values[P[i]] += pe2->values[i];
     85                this->gglobaldoflist[P[i]]=pe2->gglobaldoflist[i];
     86        }
     87
     88        /*Fset*/
     89        fsize=pe1->fsize;
     90        for(i=0;i<pe2->fsize;i++){
     91                if(P[pe2->flocaldoflist[i]] >= pe1->nrows) fsize++;
     92        }
     93        this->fsize=fsize;
     94        if(fsize){
     95                this->flocaldoflist =(int*)xmalloc(fsize*sizeof(int));
     96                this->fglobaldoflist=(int*)xmalloc(fsize*sizeof(int));
     97                for(i=0;i<pe1->fsize;i++){
     98                        this->flocaldoflist[i] =pe1->flocaldoflist[i];
     99                        this->fglobaldoflist[i]=pe1->fglobaldoflist[i];
     100                }
     101                counter=pe1->fsize;
     102                for(i=0;i<pe2->fsize;i++){
     103                        if(P[pe2->flocaldoflist[i]] >= pe1->nrows){
     104                                this->flocaldoflist[counter] =P[pe2->flocaldoflist[i]];
     105                                this->fglobaldoflist[counter]=pe2->fglobaldoflist[i];
     106                        }
     107                }
     108        }
     109        else{
     110                this->flocaldoflist=NULL;
     111                this->fglobaldoflist=NULL;
     112        }
     113
     114        /*clean-up*/
     115        xfree((void**)&P);
     116}
     117/*}}}*/
     118/*FUNCTION ElementVector::ElementVector(ElementVector* pe1, ElementVector* pe2,ElementVector* pe3){{{1*/
     119ElementVector::ElementVector(ElementVector* pe1, ElementVector* pe2,ElementVector* pe3){
     120
     121        /*Concatenate all matrices*/
     122        ElementVector* pe12 =new ElementVector(pe1,pe2);
     123        ElementVector* pe123=new ElementVector(pe12,pe3);
     124
     125        /*Initialize current object with this matrix*/
     126        this->Init(pe123);
     127
     128        /*clean-up*/
     129        delete pe12;
     130        delete pe123;
     131}
     132/*}}}*/
    33133/*FUNCTION ElementVector::ElementVector(int gsize,int* gglobaldoflist){{{1*/
    34134ElementVector::ElementVector(int gsize,int* in_gglobaldoflist){
    … …  
    49149        }
    50150        /*not needed: */
     151        this->fsize=0;
    51152        this->flocaldoflist=NULL;
    52153        this->fglobaldoflist=NULL;
    … …  
    130231}
    131232/*}}}*/
     233/*FUNCTION ElementVector::Echo{{{1*/
     234void ElementVector::Echo(void){
     235
     236        int i,j;
     237        printf("Element Vector echo: \n");
     238        printf("   nrows: %i\n",nrows);
     239        printf("   pf: %s\n",pf?"true":"false");
     240
     241        printf("   values: \n");
     242        for(i=0;i<nrows;i++){
     243                printf("      %i: %10g\n",i,values[i]);
     244        }
     245
     246        printf("   gglobaldoflist (%p): ",gglobaldoflist);
     247        if(gglobaldoflist) for(i=0;i<nrows;i++)printf("%i ",gglobaldoflist[i]); printf("\n");
     248
     249        printf("   fsize: %i\n",fsize);
     250        printf("   flocaldoflist (%p): ",flocaldoflist);
     251        if(flocaldoflist) for(i=0;i<fsize;i++)printf("%i ",flocaldoflist[i]); printf("\n");
     252        printf("   fglobaldoflist (%p): ",fglobaldoflist);
     253        if(fglobaldoflist)for(i=0;i<fsize;i++)printf("%i ",fglobaldoflist[i]); printf("\n");
     254}
     255/*}}}*/
     256/*FUNCTION ElementVector::Init{{{1*/
     257void ElementVector::Init(ElementVector* pe){
     258
     259        ISSMASSERT(pe);
     260
     261        this->nrows =pe->nrows;
     262        this->pf    =pe->pf;
     263
     264        this->values=(double*)xmalloc(this->nrows*sizeof(double));
     265        memcpy(this->values,pe->values,this->nrows*sizeof(double));
     266
     267        this->gglobaldoflist=(int*)xmalloc(this->nrows*sizeof(int));
     268        memcpy(this->gglobaldoflist,pe->gglobaldoflist,this->nrows*sizeof(int));
     269
     270        this->fsize=pe->fsize;
     271        if(this->fsize){
     272                this->flocaldoflist=(int*)xmalloc(this->fsize*sizeof(int));
     273                memcpy(this->flocaldoflist,pe->flocaldoflist,this->fsize*sizeof(int));
     274                this->fglobaldoflist=(int*)xmalloc(this->fsize*sizeof(int));
     275                memcpy(this->fglobaldoflist,pe->fglobaldoflist,this->fsize*sizeof(int));
     276        }
     277        else{
     278                this->flocaldoflist=NULL;
     279                this->fglobaldoflist=NULL;
     280        }
     281}
     282/*}}}*/

  • TabularUnified issm/trunk/src/c/objects/Numerics/ElementVector.h

    r5827 r5926  
    3333                /*ElementVector constructors, destructors {{{1*/
    3434                ElementVector();
     35                ElementVector(ElementVector* pe1,ElementVector* pe2);
     36                ElementVector(ElementVector* pe1,ElementVector* pe2,ElementVector* pe3);
    3537                ElementVector(int gsize,int* gglobaldoflist);
    3638                ElementVector(int gsize,int* flocaldoflist,int* fglobaldoflist,int fsize);
    … …  
    4042                void AddValues(double* pe_gg);
    4143                void AddToGlobal(Vec pg, Vec pf);
     44                void Echo(void);
     45                void Init(ElementVector* pe);
    4246                /*}}}*/
    4347};
Note: See TracChangeset for help on using the changeset viewer.