Changeset 16926


Ignore:
Timestamp:
11/25/13 11:28:46 (11 years ago)
Author:
seroussi
Message:

CHG: keep working on coupling PVectors

Location:
issm/trunk-jpl/src/c
Files:
8 edited

Legend:

Unmodified
Added
Removed

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

    r16922 r16926  
    851851                case SSAHOApproximationEnum:
    852852                        return CreatePVectorSSAHO(element);
     853                case HOFSApproximationEnum:
     854                        return CreatePVectorHOFS(element);
    853855                case NoneApproximationEnum:
    854856                        return NULL;
    … …  
    31013103        int         i,j;
    31023104        IssmDouble  Jdet2d,alpha2;
    3103         IssmDouble *xyz_list_tria=NULL;
    3104         IssmDouble* L           = xNewZeroInit<IssmDouble>(2*numdof);
    3105         IssmDouble  DL[2][2]     ={{ 0,0 },{0,0}}; //for basal drag
     3105        IssmDouble *xyz_list_tria = NULL;
     3106        IssmDouble* L             = xNewZeroInit<IssmDouble>(2*numdof);
     3107        IssmDouble  DL[2][2]      = {{ 0,0 },{0,0}}; //for basal drag
    31063108        IssmDouble  DL_scalar;
    3107         IssmDouble* Ke_gg       = xNewZeroInit<IssmDouble>(numdof*numdof);
    3108         Node       *node_list[2*numnodes];
    3109         int*        cs_list= xNew<int>(2*numnodes);
     3109        IssmDouble* Ke_gg         = xNewZeroInit<IssmDouble>(numdof*numdof);
     3110        Node      **node_list     = xNew<Node*>(2*numnodes);
     3111        int*        cs_list       = xNew<int>(2*numnodes);
    31103112
    31113113        /*Initialize Element matrix and return if necessary*/
    … …  
    31633165        delete friction;
    31643166        xDelete<int>(cs_list);
     3167        xDelete<Node*>(node_list);
    31653168        xDelete<IssmDouble>(xyz_list_tria);
    31663169        xDelete<IssmDouble>(Ke_gg);
    … …  
    31843187        IssmDouble  D_scalar;
    31853188        IssmDouble* Ke_gg          = xNewZeroInit<IssmDouble>(numdofp*numdofm);
    3186         Node        *node_list[2*numnodes];
     3189        Node       **node_list     = xNew<Node*>(2*numnodes);
    31873190        int*         cs_list= xNew<int>(2*numnodes);
    31883191
    … …  
    32493252        xDelete<IssmDouble>(Ke_gg);
    32503253        xDelete<IssmDouble>(xyz_list);
     3254        xDelete<Node*>(node_list);
    32513255        xDelete<int>(cs_list);
    32523256        return Ke;
    … …  
    32653269        return pe;
    32663270}/*}}}*/
     3271ElementVector* StressbalanceAnalysis::CreatePVectorHOFS(Element* element){/*{{{*/
     3272
     3273        /*compute all load vectors for this element*/
     3274        int init = element->FiniteElement();
     3275        element->SetTemporaryElementType(P1Enum);
     3276        ElementVector* pe1=CreatePVectorHO(element);
     3277        element->SetTemporaryElementType(init);
     3278        ElementVector* pe2=CreatePVectorFS(element);
     3279        int indices[3]={18,19,20};
     3280        element->SetTemporaryElementType(MINIcondensedEnum);
     3281        ElementMatrix* Ke = CreateKMatrixFS(element);
     3282        element->SetTemporaryElementType(init);
     3283        pe2->StaticCondensation(Ke,3,&indices[0]);
     3284        delete Ke;
     3285        ElementVector* pe3=CreatePVectorCouplingHOFS(element);
     3286        ElementVector* pe =new ElementVector(pe1,pe2,pe3);
     3287
     3288        /*clean-up and return*/
     3289        delete pe1;
     3290        delete pe2;
     3291        delete pe3;
     3292        return pe;
     3293}/*}}}*/
     3294ElementVector* StressbalanceAnalysis::CreatePVectorCouplingHOFS(Element* element){/*{{{*/
     3295
     3296        /*compute all load vectors for this element*/
     3297        ElementVector* pe1=CreatePVectorCouplingHOFSViscous(element);
     3298        ElementVector* pe2=CreatePVectorCouplingHOFSFriction(element);
     3299        ElementVector* pe =new ElementVector(pe1,pe2);
     3300
     3301        /*clean-up and return*/
     3302        delete pe1;
     3303        delete pe2;
     3304        return pe;
     3305}/*}}}*/
     3306ElementVector* StressbalanceAnalysis::CreatePVectorCouplingHOFSFriction(Element* element){/*{{{*/
     3307
     3308        /*Intermediaries*/
     3309        int         i,j,approximation;
     3310        IssmDouble  Jdet,Jdet2d,FSreconditioning;
     3311        IssmDouble      bed_normal[3];
     3312        IssmDouble  viscosity, w, alpha2_gauss;
     3313        IssmDouble  dw[3];
     3314        IssmDouble      *xyz_list_tria = NULL;
     3315        IssmDouble  *xyz_list      = NULL;
     3316        IssmDouble  basis[6]; //for the six nodes of the penta
     3317
     3318        /*Initialize Element vector and return if necessary*/
     3319        if(!element->IsOnBed() || element->IsFloating()) return NULL;
     3320        element->GetInputValue(&approximation,ApproximationEnum);
     3321        if(approximation!=HOFSApproximationEnum) return NULL;
     3322
     3323        int vnumnodes = element->GetNumberOfNodesVelocity();
     3324        int pnumnodes = element->GetNumberOfNodesPressure();
     3325        int numnodes  = vnumnodes+pnumnodes;
     3326
     3327        /*Prepare coordinate system list*/
     3328        int*   cs_list   = xNew<int>(vnumnodes+pnumnodes);
     3329        Node **node_list = xNew<Node*>(vnumnodes+pnumnodes);
     3330        for(i=0;i<vnumnodes;i++){
     3331                cs_list[i]           = XYZEnum;
     3332                node_list[i]           = element->GetNode(i);
     3333        }
     3334        for(i=0;i<pnumnodes;i++){
     3335                cs_list[vnumnodes+i] = PressureEnum;
     3336                node_list[vnumnodes+i] = element->GetNode(vnumnodes+i);
     3337        }
     3338
     3339        ElementVector* pe=element->NewElementVector(FSvelocityEnum);
     3340
     3341        /*Retrieve all inputs and parameters*/
     3342        element->GetVerticesCoordinates(&xyz_list);
     3343        element->GetVerticesCoordinatesBase(&xyz_list_tria);
     3344        element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
     3345        Input* vx_input=  element->GetInput(VxEnum);   _assert_(vx_input);
     3346        Input* vy_input=  element->GetInput(VyEnum);   _assert_(vy_input);
     3347        Input* vz_input=  element->GetInput(VzEnum);   _assert_(vz_input);
     3348        Input* vzHO_input=element->GetInput(VzHOEnum); _assert_(vzHO_input);
     3349
     3350        /*build friction object, used later on: */
     3351        Friction* friction=new Friction(element,3);
     3352
     3353        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
     3354        Gauss* gauss=element->NewGaussBase(2);
     3355        for(int ig=gauss->begin();ig<gauss->end();ig++){
     3356
     3357                gauss->GaussPoint(ig);
     3358
     3359                element->JacobianDeterminantBase(&Jdet2d,xyz_list_tria,gauss);
     3360                element->NodalFunctionsP1(basis, gauss);
     3361
     3362                vzHO_input->GetInputValue(&w, gauss);
     3363                vzHO_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss);
     3364
     3365                element->NormalBase(&bed_normal[0],xyz_list_tria);
     3366                element->ViscosityFS(&viscosity,xyz_list,gauss,vx_input,vy_input,vz_input);
     3367                friction->GetAlpha2(&alpha2_gauss, gauss,vx_input,vy_input,vz_input);
     3368
     3369                for(i=0;i<3;i++){
     3370                        pe->values[i*3+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i];
     3371                        pe->values[i*3+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i];
     3372                        pe->values[i*3+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*basis[i];
     3373                }
     3374        }
     3375
     3376        /*Transform coordinate system*/
     3377        element->TransformLoadVectorCoord(pe,node_list,vnumnodes+pnumnodes,cs_list);
     3378
     3379        pe->Echo();
     3380        element->CreatePVector();
     3381        /*Clean up and return*/
     3382        xDelete<int>(cs_list);
     3383        xDelete<Node*>(node_list);
     3384        xDelete<IssmDouble>(xyz_list);
     3385        xDelete<IssmDouble>(xyz_list_tria);
     3386        delete gauss;
     3387        delete friction;
     3388        return pe;
     3389}/*}}}*/
     3390ElementVector* StressbalanceAnalysis::CreatePVectorCouplingHOFSViscous(Element* element){/*{{{*/
     3391
     3392        /*Intermediaries */
     3393        int         i,approximation;
     3394        IssmDouble  viscosity,Jdet,FSreconditioning;
     3395        IssmDouble  dw[3];
     3396        IssmDouble  *xyz_list = NULL;
     3397        IssmDouble  basis[6]; //for the six nodes of the penta
     3398        IssmDouble  dbasis[3][6]; //for the six nodes of the penta
     3399
     3400        /*Initialize Element vector and return if necessary*/
     3401        element->GetInputValue(&approximation,ApproximationEnum);
     3402        if(approximation!=HOFSApproximationEnum) return NULL;
     3403        int   vnumnodes = element->GetNumberOfNodesVelocity();
     3404        int   pnumnodes = element->GetNumberOfNodesPressure();
     3405
     3406        /*Prepare coordinate system list*/
     3407        int*   cs_list   = xNew<int>(vnumnodes+pnumnodes);
     3408        Node **node_list = xNew<Node*>(vnumnodes+pnumnodes);
     3409        for(i=0;i<vnumnodes;i++){
     3410                cs_list[i]             = XYZEnum;
     3411                node_list[i]           = element->GetNode(i);
     3412        }
     3413        for(i=0;i<pnumnodes;i++){
     3414                cs_list[vnumnodes+i]   = PressureEnum;
     3415                node_list[vnumnodes+i] = element->GetNode(vnumnodes+i);
     3416        }
     3417        ElementVector* pe = element->NewElementVector(FSvelocityEnum);
     3418
     3419        /*Retrieve all inputs and parameters*/
     3420        element->GetVerticesCoordinates(&xyz_list);
     3421        element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
     3422        Input* vx_input   =element->GetInput(VxEnum);   _assert_(vx_input);
     3423        Input* vy_input   =element->GetInput(VyEnum);   _assert_(vy_input);
     3424        Input* vz_input   =element->GetInput(VzEnum);   _assert_(vz_input);
     3425        Input* vzHO_input=element->GetInput(VzHOEnum);  _assert_(vzHO_input);
     3426
     3427        /* Start  looping on the number of gaussian points: */
     3428        Gauss* gauss=element->NewGauss(5);
     3429        for(int ig=gauss->begin();ig<gauss->end();ig++){
     3430
     3431                gauss->GaussPoint(ig);
     3432
     3433                element->JacobianDeterminant(&Jdet, xyz_list,gauss);
     3434                element->NodalFunctionsP1(&basis[0],gauss);
     3435                element->NodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
     3436               
     3437                element->ViscosityFS(&viscosity,xyz_list,gauss,vx_input,vy_input,vz_input);
     3438                vzHO_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss);
     3439
     3440                for(i=0;i<vnumnodes;i++){
     3441                        pe->values[i*3+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i];
     3442                        pe->values[i*3+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i];
     3443                        pe->values[i*3+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]);
     3444                        pe->values[3*vnumnodes+i]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i];
     3445                }
     3446        }
     3447
     3448        /*Transform coordinate system*/
     3449        element->TransformLoadVectorCoord(pe,node_list,vnumnodes+pnumnodes,cs_list);
     3450
     3451        /*Clean up and return*/
     3452        xDelete<int>(cs_list);
     3453        xDelete<Node*>(node_list);
     3454        xDelete<IssmDouble>(xyz_list);
     3455        delete gauss;
     3456        return pe;
     3457}/*}}}*/
     3458
    32673459void StressbalanceAnalysis::GetBSSAHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    32683460        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF2.

  • TabularUnified issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h

    r16907 r16926  
    7676                ElementMatrix* CreateKMatrixCouplingSSAHOViscous(Element* element);
    7777                ElementVector* CreatePVectorSSAHO(Element* element);
     78                ElementVector* CreatePVectorHOFS(Element* element);
     79                ElementVector* CreatePVectorCouplingHOFS(Element* element);
     80                ElementVector* CreatePVectorCouplingHOFSFriction(Element* element);
     81                ElementVector* CreatePVectorCouplingHOFSViscous(Element* element);
    7882                void GetBSSAHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    7983                void InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element);

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

    r16922 r16926  
    314314        ::TransformSolutionCoord(values,this->nodes,numnodes,transformenum_list);
    315315}/*}}}*/
     316void Element::TransformLoadVectorCoord(ElementVector* Ke,Node** nodes_list,int numnodes,int* transformenum_list){/*{{{*/
     317        ::TransformLoadVectorCoord(Ke,nodes_list,numnodes,transformenum_list);
     318}/*}}}*/
    316319void Element::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum){/*{{{*/
    317320        ::TransformStiffnessMatrixCoord(Ke,this->nodes,this->GetNumberOfNodes(),transformenum);

  • TabularUnified issm/trunk-jpl/src/c/classes/Elements/Element.h

    r16922 r16926  
    7777                void       TransformLoadVectorCoord(ElementVector* pe,int numnodes,int transformenum){_error_("not implemented yet");};      /*Tiling only*/
    7878                void       TransformLoadVectorCoord(ElementVector* pe,int numnodes,int* transformenum_list){_error_("not implemented yet");};      /*Tiling only*/
     79                void       TransformLoadVectorCoord(ElementVector* pe,Node** nodes_list,int numnodes,int* transformenum_list);
    7980                void       TransformSolutionCoord(IssmDouble* values,int transformenum);
    8081                void       TransformSolutionCoord(IssmDouble* values,int* transformenum_list);
    … …  
    106107                virtual IssmDouble MinEdgeLength(IssmDouble* xyz_list)=0;
    107108                virtual void   NodalFunctions(IssmDouble* basis,Gauss* gauss)=0;
     109                virtual void   NodalFunctionsP1(IssmDouble* basis,Gauss* gauss)=0;
    108110                virtual void   NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss)=0;
    109111                virtual void   NodalFunctionsPressure(IssmDouble* basis, Gauss* gauss)=0;
    110112                virtual void   NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss)=0;
     113                virtual void   NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss)=0;
    111114                virtual void   NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss)=0;
    112115                virtual void   NormalSection(IssmDouble* normal,IssmDouble* xyz_list)=0;
    … …  
    193196                virtual void   ReduceMatrices(ElementMatrix* Ke,ElementVector* pe)=0;
    194197                virtual void   ResetCoordinateSystem()=0;
     198                virtual void   SetTemporaryElementType(int element_type_in)=0;
    195199                virtual IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa)=0;
    196200                virtual void   ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss)=0;

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

    r16921 r16926  
    25062506}
    25072507/*}}}*/
     2508/*FUNCTION Penta::NodalFunctionsP1{{{*/
     2509void Penta::NodalFunctionsP1(IssmDouble* basis, Gauss* gauss){
     2510
     2511        _assert_(gauss->Enum()==GaussPentaEnum);
     2512        this->GetNodalFunctionsP1(basis,(GaussPenta*)gauss);
     2513
     2514}
     2515/*}}}*/
    25082516/*FUNCTION Penta::NodalFunctionsDerivatives{{{*/
    25092517void Penta::NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){
    … …  
    25112519        _assert_(gauss->Enum()==GaussPentaEnum);
    25122520        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss);
     2521
     2522}
     2523/*}}}*/
     2524/*FUNCTION Penta::NodalFunctionsP1Derivatives{{{*/
     2525void Penta::NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){
     2526
     2527        _assert_(gauss->Enum()==GaussPentaEnum);
     2528        this->GetNodalFunctionsP1Derivatives(dbasis,xyz_list,(GaussPenta*)gauss);
    25132529
    25142530}
    … …  
    28472863        if (this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
    28482864        else this->nodes=NULL;
     2865}
     2866/*}}}*/
     2867/*FUNCTION Penta::SetTemporaryElementType{{{*/
     2868void Penta::SetTemporaryElementType(int element_type_in){
     2869        this->element_type=element_type_in;
    28492870}
    28502871/*}}}*/
    … …  
    82348255        ::TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list);
    82358256
     8257        pe->Echo();
     8258        _error_("stop");
    82368259        /*Clean up and return*/
    82378260        xDelete<int>(cs_list);
    … …  
    83968419        /*compute all load vectors for this element*/
    83978420        int init = this->element_type;
     8421        printf("init init %i\n",init);
    83988422        this->element_type=P1Enum;
    83998423        ElementVector* pe1=CreatePVectorStressbalanceHO();

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

    r16916 r16926  
    115115                void   ResetCoordinateSystem(void);
    116116                void   SmbGradients();
     117                void   SetTemporaryElementType(int element_type_in);
    117118           Element* SpawnBasalElement(void);
    118119                IssmDouble SurfaceArea(void);
    … …  
    257258                ElementMatrix* NewElementMatrix(int approximation_enum);
    258259                void           NodalFunctions(IssmDouble* basis,Gauss* gauss);
     260                void           NodalFunctionsP1(IssmDouble* basis,Gauss* gauss);
    259261                void           NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
     262                void           NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
    260263                void           NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
    261264                void           NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss);

  • TabularUnified issm/trunk-jpl/src/c/classes/Elements/Seg.h

    r16916 r16926  
    108108                IssmDouble  MinEdgeLength(IssmDouble* xyz_list){_error_("not implemented yet");};
    109109                void        NodalFunctions(IssmDouble* basis,Gauss* gauss);
     110                void        NodalFunctionsP1(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");};
    110111                void        NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");};
    111112                void        NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");};
    112113                void        NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
     114                void        NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    113115                void        NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    114116                bool        NoIceInElement(){_error_("not implemented yet");};
    … …  
    181183                void        ResetCoordinateSystem(void){_error_("not implemented yet");};
    182184                void        ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){_error_("not implemented yet");};
     185                void        SetTemporaryElementType(int element_type_in){_error_("not implemented yet");};
    183186                void          SmbGradients(){_error_("not implemented yet");};
    184187                IssmDouble  SurfaceArea(void){_error_("not implemented yet");};

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

    r16916 r16926  
    288288                ElementMatrix* NewElementMatrix(int approximation_enum);
    289289                void           NodalFunctions(IssmDouble* basis,Gauss* gauss);
     290                void           NodalFunctionsP1(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");};
    290291                void           NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
     292                void           NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    291293                void           NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
    292294                void           NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss);
    293295                void           NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss);
    294296                void             SetClone(int* minranks);
     297                void           SetTemporaryElementType(int element_type_in){_error_("not implemented yet");};
    295298                Seg*             SpawnSeg(int index1,int index2);
    296299                IssmDouble     StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");};
Note: See TracChangeset for help on using the changeset viewer.