Changeset 25610


Ignore:
Timestamp:
09/29/20 11:24:39 (5 years ago)
Author:
Mathieu Morlighem
Message:

CHG: working on MLHO

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

Legend:

Unmodified
Added
Removed

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

    r25554 r25610  
    3131        IssmDouble rho_ice;
    3232        IssmDouble FSreconditioning;
    33         bool       isSIA,isSSA,isL1L2,isHO,isFS,iscoupling;
     33        bool       isSIA,isSSA,isL1L2,isMLHO,isHO,isFS,iscoupling;
    3434        bool       spcpresent = false;
    3535        int        Mx,Nx;
    … …  
    5959        iomodel->FindConstant(&isSSA,"md.flowequation.isSSA");
    6060        iomodel->FindConstant(&isL1L2,"md.flowequation.isL1L2");
     61        iomodel->FindConstant(&isMLHO,"md.flowequation.isMLHO");
    6162        iomodel->FindConstant(&isHO,"md.flowequation.isHO");
    6263        iomodel->FindConstant(&isFS,"md.flowequation.isFS");
    6364
    64         /*Now, is the flag macayaealHO on? otherwise, do nothing: */
    65         if(!isSSA && !isHO && !isFS && !isL1L2) return;
     65        /*Is this model only SIA??*/
     66        if(!isSSA && !isHO && !isFS && !isL1L2 && !isMLHO) return;
    6667
    6768        /*Do we have coupling*/
    68         if((isSIA?1.:0.) + (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
     69        if((isSIA?1.:0.) + (isSSA?1.:0.) + (isL1L2?1.:0.) + (isMLHO?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
    6970         iscoupling = true;
    7071        else
    … …  
    7778                if(isSSA)       iomodel->FindConstant(&finiteelement,"md.flowequation.fe_SSA");
    7879                else if(isL1L2) finiteelement = P1Enum;
     80                else if(isMLHO) finiteelement = P1Enum;
    7981                else if(isHO)   iomodel->FindConstant(&finiteelement,"md.flowequation.fe_HO");
    8082                else if(isFS){  iomodel->FindConstant(&finiteelement,"md.flowequation.fe_FS");
    … …  
    452454        int         count;
    453455        int         penpair_ids[2];
    454         bool        isSSA,isL1L2,isHO,isFS;
     456        bool        isSSA,isL1L2,isMLHO,isHO,isFS;
    455457        int         numpenalties,numrifts,numriftsegments;
    456458        IssmDouble *riftinfo       = NULL;
    … …  
    460462        /*Fetch parameters: */
    461463        iomodel->FindConstant(&isL1L2,"md.flowequation.isL1L2");
     464        iomodel->FindConstant(&isMLHO,"md.flowequation.isMLHO");
    462465        iomodel->FindConstant(&isFS,"md.flowequation.isFS");
    463466        iomodel->FindConstant(&isSSA,"md.flowequation.isSSA");
    … …  
    465468        iomodel->FindConstant(&numrifts,"md.rifts.numrifts");
    466469
    467         /*Now, is the flag macayaealHO on? otherwise, do nothing: */
    468         if(!isSSA && !isHO && !isFS && !isL1L2) return;
     470        /*Is this SIA only?*/
     471        if(!isSSA && !isHO && !isFS && !isL1L2 && !isMLHO) return;
    469472
    470473        /*Initialize counter: */
    … …  
    511514
    512515        /*Intermediary*/
    513         bool isSSA,isL1L2,isHO,isFS,iscoupling;
     516        bool isSSA,isL1L2,isMLHO,isHO,isFS,iscoupling;
    514517        int  finiteelement=-1,approximation=-1;
    515518
    … …  
    517520        iomodel->FindConstant(&isSSA,"md.flowequation.isSSA");
    518521        iomodel->FindConstant(&isL1L2,"md.flowequation.isL1L2");
     522        iomodel->FindConstant(&isMLHO,"md.flowequation.isMLHO");
    519523        iomodel->FindConstant(&isHO,"md.flowequation.isHO");
    520524        iomodel->FindConstant(&isFS,"md.flowequation.isFS");
    521525
    522526        /*Now, check that we have non SIA elements */
    523         if(!isSSA & !isL1L2 & !isHO & !isFS) return;
     527        if(!isSSA && !isL1L2 && !isMLHO && !isHO && !isFS) return;
    524528
    525529        /*Do we have coupling*/
    526         if( (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
     530        if( (isSSA?1.:0.) + (isL1L2?1.:0.) + (isMLHO?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
    527531         iscoupling = true;
    528532        else
    … …  
    539543                else if(isL1L2){
    540544                        approximation = L1L2ApproximationEnum;
     545                        finiteelement = P1Enum;
     546                }
     547                else if(isMLHO){
     548                        approximation = MLHOApproximationEnum;
    541549                        finiteelement = P1Enum;
    542550                }
    … …  
    613621                         }
    614622                         break;
    615                 case L1L2ApproximationEnum: numdofs =2; break;
     623                case L1L2ApproximationEnum: numdofs = 2; break;
     624                case MLHOApproximationEnum: numdofs = 4; break;
    616625                case HOApproximationEnum:
    617626                         switch(domaintype){
    … …  
    678687        int    FrictionCoupling;
    679688        int*   finiteelement_list=NULL;
    680         bool   isSSA,isL1L2,isHO,isFS,iscoupling;
     689        bool   isSSA,isL1L2,isMLHO,isHO,isFS,iscoupling;
    681690        bool   control_analysis;
    682691        bool   dakota_analysis;
    … …  
    686695        iomodel->FindConstant(&isSSA,"md.flowequation.isSSA");
    687696        iomodel->FindConstant(&isL1L2,"md.flowequation.isL1L2");
     697        iomodel->FindConstant(&isMLHO,"md.flowequation.isMLHO");
    688698        iomodel->FindConstant(&isHO,"md.flowequation.isHO");
    689699        iomodel->FindConstant(&isFS,"md.flowequation.isFS");
    … …  
    695705
    696706        /*return if no processing required*/
    697         if(!isSSA & !isL1L2 & !isHO & !isFS) return;
     707        if(!isSSA && !isL1L2 && !isMLHO && !isHO && !isFS) return;
    698708
    699709        /*Fetch data needed and allocate vectors: */
    … …  
    702712
    703713        /*Do we have coupling*/
    704         if( (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
     714        if( (isSSA?1.:0.) + (isL1L2?1.:0.) + (isMLHO?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
    705715         iscoupling = true;
    706716        else
    … …  
    711721                if(isSSA)       iomodel->FindConstant(&finiteelement,"md.flowequation.fe_SSA");
    712722                else if(isL1L2) finiteelement = P1Enum;
     723                else if(isMLHO) finiteelement = P1Enum;
    713724                else if(isHO)   iomodel->FindConstant(&finiteelement,"md.flowequation.fe_HO");
    714725                else if(isFS)   iomodel->FindConstant(&finiteelement,"md.flowequation.fe_FS");
    … …  
    926937        parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.isSSA",FlowequationIsSSAEnum));
    927938        parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.isL1L2",FlowequationIsL1L2Enum));
     939        parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.isMLHO",FlowequationIsMLHOEnum));
    928940        parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.isHO",FlowequationIsHOEnum));
    929941        parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.isFS",FlowequationIsFSEnum));
    … …  
    10231035
    10241036        /*Intermediaries*/
    1025         bool isSSA,isL1L2,isHO,isFS;
     1037        bool isSSA,isL1L2,isMLHO,isHO,isFS;
    10261038        bool conserve_loads = true;
    10271039        int  newton,domaintype,fe_FS;
    … …  
    10301042        femmodel->parameters->FindParam(&isSSA,FlowequationIsSSAEnum);
    10311043        femmodel->parameters->FindParam(&isL1L2,FlowequationIsL1L2Enum);
     1044        femmodel->parameters->FindParam(&isMLHO,FlowequationIsMLHOEnum);
    10321045        femmodel->parameters->FindParam(&isHO,FlowequationIsHOEnum);
    10331046        femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
    … …  
    10361049        femmodel->parameters->FindParam(&newton,StressbalanceIsnewtonEnum);
    10371050
    1038         if(isFS && !(isSSA || isHO || isL1L2)){
     1051        if(isFS && !(isSSA || isHO || isL1L2 || isMLHO)){
    10391052                if(VerboseSolution()) _printf0_("   computing velocities\n");
    10401053
    … …  
    10611074                 solutionsequence_nonlinear(femmodel,conserve_loads);
    10621075        }
    1063         else if(!isFS && (isSSA || isHO || isL1L2)){
     1076        else if(!isFS && (isSSA || isHO || isL1L2 || isMLHO)){
    10641077                if(VerboseSolution()) _printf0_("   computing velocities\n");
    10651078
    … …  
    10771090                }
    10781091        }
    1079         else if ((isSSA || isL1L2 || isHO) && isFS){
     1092        else if ((isSSA || isL1L2 || isMLHO || isHO) && isFS){
    10801093                if(VerboseSolution()) _printf0_("   computing coupling between lower order models and FS\n");
    10811094                solutionsequence_FScoupling_nonlinear(femmodel,conserve_loads);
    … …  
    11261139                case L1L2ApproximationEnum:
    11271140                        return CreateKMatrixL1L2(element);
     1141                case MLHOApproximationEnum:
     1142                        return CreateKMatrixMLHO(element);
    11281143                case HOApproximationEnum:
    11291144                        return CreateKMatrixHO(element);
    … …  
    11531168                case L1L2ApproximationEnum:
    11541169                        return CreatePVectorL1L2(element);
     1170                case MLHOApproximationEnum:
     1171                        return CreatePVectorMLHO(element);
    11551172                case HOApproximationEnum:
    11561173                        return CreatePVectorHO(element);
    … …  
    11771194                        GetSolutionFromInputsFS(solution,element);
    11781195                        return;
    1179                 case SSAApproximationEnum: case HOApproximationEnum: case SIAApproximationEnum:
    1180                         GetSolutionFromInputsHoriz(solution,element);
    1181                         return;
    1182                 case L1L2ApproximationEnum:
     1196                case SSAApproximationEnum: case HOApproximationEnum: case L1L2ApproximationEnum: case MLHOApproximationEnum: case SIAApproximationEnum:
    11831197                        GetSolutionFromInputsHoriz(solution,element);
    11841198                        return;
    … …  
    12611275                case L1L2ApproximationEnum:
    12621276                        InputUpdateFromSolutionL1L2(solution,element);
     1277                        return;
     1278                case MLHOApproximationEnum:
     1279                        InputUpdateFromSolutionMLHO(solution,element);
    12631280                        return;
    12641281                case SSAHOApproximationEnum:
    … …  
    23882405        /*Fetch number of nodes and dof for this finite element*/
    23892406        int numnodes = basalelement->GetNumberOfNodes();
    2390         int numdof   = numnodes*2;
    23912407
    23922408        /*Initialize Element matrix and vectors*/
    23932409        ElementMatrix* Ke     = basalelement->NewElementMatrix(L1L2ApproximationEnum);
    2394         IssmDouble*    B      = xNew<IssmDouble>(3*numdof);
    2395         IssmDouble*    Bprime = xNew<IssmDouble>(3*numdof);
    2396         IssmDouble*    D      = xNewZeroInit<IssmDouble>(3*3);
     2410        IssmDouble*    dbasis = xNew<IssmDouble>(3*numnodes);
    23972411
    23982412        /*Retrieve all inputs and parameters*/
    … …  
    24092423
    24102424                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    2411                 this->GetBSSA(B,basalelement,2,xyz_list,gauss_base);
    2412                 this->GetBSSAprime(Bprime,basalelement,2,xyz_list,gauss_base);
     2425                basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss_base);
    24132426
    24142427                element->material->ViscosityL1L2(&viscosity,xyz_list,gauss,vx_input,vy_input,surface_input);
    24152428
    2416                 for(int i=0;i<3;i++) D[i*3+i]=2*viscosity*gauss->weight*Jdet;
    2417 
    2418                 TripleMultiply(B,3,numdof,1,
    2419                                         D,3,3,0,
    2420                                         Bprime,3,numdof,0,
    2421                                         &Ke->values[0],1);
     2429                for(int i=0;i<numnodes;i++){
     2430                        for(int j=0;j<numnodes;j++){
     2431                                Ke->values[2*i*2*numnodes+2*j] += gauss->weight*Jdet*viscosity*(
     2432                                                        4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
     2433                                                        );
     2434                                Ke->values[2*i*2*numnodes+2*j+1] += gauss->weight*Jdet*viscosity*(
     2435                                                        2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i]
     2436                                                        );
     2437                                Ke->values[(2*i+1)*2*numnodes+2*j] += gauss->weight*Jdet*viscosity*(
     2438                                                        2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i]
     2439                                                        );
     2440                                Ke->values[(2*i+1)*2*numnodes+2*j+1] += gauss->weight*Jdet*viscosity*(
     2441                                                        dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + 4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
     2442                                                        );
     2443                        }
     2444                }
    24222445        }
    24232446
    … …  
    24302453        basalelement->DeleteMaterials(); delete basalelement;
    24312454        xDelete<IssmDouble>(xyz_list);
    2432         xDelete<IssmDouble>(D);
    2433         xDelete<IssmDouble>(Bprime);
    2434         xDelete<IssmDouble>(B);
     2455        xDelete<IssmDouble>(dbasis);
    24352456        return Ke;
    24362457}/*}}}*/
    … …  
    25732594}/*}}}*/
    25742595void           StressbalanceAnalysis::InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element){/*{{{*/
     2596
     2597        int         i,dim,domaintype;
     2598        IssmDouble  rho_ice,g;
     2599        int*        doflist=NULL;
     2600        IssmDouble* xyz_list=NULL;
     2601        Element*    basalelement=NULL;
     2602
     2603        /*Deal with pressure first*/
     2604        int numvertices = element->GetNumberOfVertices();
     2605        IssmDouble* pressure  = xNew<IssmDouble>(numvertices);
     2606        IssmDouble* thickness = xNew<IssmDouble>(numvertices);
     2607        IssmDouble* surface   = xNew<IssmDouble>(numvertices);
     2608
     2609        element->FindParam(&dim,DomainDimensionEnum);
     2610        element->FindParam(&domaintype,DomainTypeEnum);
     2611        rho_ice =element->FindParam(MaterialsRhoIceEnum);
     2612        g       =element->FindParam(ConstantsGEnum);
     2613        if(dim==2){
     2614                element->GetInputListOnVertices(thickness,ThicknessEnum);
     2615                for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i];
     2616        }
     2617        else{
     2618                element->GetVerticesCoordinates(&xyz_list);
     2619                element->GetInputListOnVertices(surface,SurfaceEnum);
     2620                for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
     2621        }
     2622        element->AddInput(PressureEnum,pressure,P1Enum);
     2623        xDelete<IssmDouble>(pressure);
     2624        xDelete<IssmDouble>(thickness);
     2625        xDelete<IssmDouble>(surface);
     2626
     2627        /*Get basal element*/
     2628        switch(domaintype){
     2629                case Domain2DhorizontalEnum:
     2630                        basalelement = element;
     2631                        break;
     2632                case Domain3DEnum:
     2633                        if(!element->IsOnBase()){xDelete<IssmDouble>(xyz_list); return;}
     2634                        basalelement=element->SpawnBasalElement();
     2635                        break;
     2636                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     2637        }
     2638
     2639        /*Fetch number of nodes and dof for this finite element*/
     2640        int numnodes = basalelement->GetNumberOfNodes();
     2641        int numdof   = numnodes*2;
     2642
     2643        /*Fetch dof list and allocate solution vectors*/
     2644        basalelement->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
     2645        IssmDouble* values    = xNew<IssmDouble>(numdof);
     2646        IssmDouble* vx        = xNew<IssmDouble>(numnodes);
     2647        IssmDouble* vy        = xNew<IssmDouble>(numnodes);
     2648        IssmDouble* vz        = xNew<IssmDouble>(numnodes);
     2649        IssmDouble* vel       = xNew<IssmDouble>(numnodes);
     2650
     2651        /*Use the dof list to index into the solution vector: */
     2652        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
     2653
     2654        /*Transform solution in Cartesian Space*/
     2655        basalelement->TransformSolutionCoord(&values[0],XYEnum);
     2656        basalelement->FindParam(&domaintype,DomainTypeEnum);
     2657
     2658        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
     2659        for(i=0;i<numnodes;i++){
     2660                vx[i]=values[i*2+0];
     2661                vy[i]=values[i*2+1];
     2662
     2663                /*Check solution*/
     2664                if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
     2665                if(xIsInf<IssmDouble>(vx[i])) _error_("Inf found in solution vector");
     2666                if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
     2667                if(xIsInf<IssmDouble>(vy[i])) _error_("Inf found in solution vector");
     2668        }
     2669
     2670        /*Get Vz and compute vel*/
     2671        basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.);
     2672        for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
     2673
     2674        /*Add vx and vy as inputs to the tria element: */
     2675        element->AddBasalInput(VxEnum,vx,element->GetElementType());
     2676        element->AddBasalInput(VyEnum,vy,element->GetElementType());
     2677        element->AddBasalInput(VelEnum,vel,element->GetElementType());
     2678
     2679        /*Free ressources:*/
     2680        xDelete<IssmDouble>(vel);
     2681        xDelete<IssmDouble>(vz);
     2682        xDelete<IssmDouble>(vy);
     2683        xDelete<IssmDouble>(vx);
     2684        xDelete<IssmDouble>(values);
     2685        xDelete<IssmDouble>(xyz_list);
     2686        xDelete<int>(doflist);
     2687        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     2688}/*}}}*/
     2689
     2690/*MLHO*/
     2691ElementMatrix* StressbalanceAnalysis::CreateKMatrixMLHO(Element* element){/*{{{*/
     2692        _error_("not implemented yet");
     2693
     2694        /* Check if ice in element */
     2695        if(!element->IsIceInElement()) return NULL;
     2696
     2697        /*compute all stiffness matrices for this element*/
     2698        ElementMatrix* Ke1=CreateKMatrixMLHOViscous(element);
     2699        ElementMatrix* Ke2=CreateKMatrixMLHOFriction(element);
     2700        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
     2701
     2702        /*clean-up and return*/
     2703        delete Ke1;
     2704        delete Ke2;
     2705        return Ke;
     2706}/*}}}*/
     2707ElementMatrix* StressbalanceAnalysis::CreateKMatrixMLHOFriction(Element* element){/*{{{*/
     2708        _error_("not implemented yet");
     2709
     2710        if(!element->IsOnBase() || element->IsFloating()) return NULL;
     2711        Element* basalelement = element->SpawnBasalElement();
     2712        ElementMatrix* Ke = CreateKMatrixSSAFriction(basalelement);
     2713
     2714        /*clean-up and return*/
     2715        basalelement->DeleteMaterials(); delete basalelement;
     2716        return Ke;
     2717}/*}}}*/
     2718ElementMatrix* StressbalanceAnalysis::CreateKMatrixMLHOViscous(Element* element){/*{{{*/
     2719        _error_("not implemented yet");
     2720
     2721        /*Intermediaries*/
     2722        IssmDouble  viscosity,Jdet;
     2723        IssmDouble *xyz_list = NULL;
     2724
     2725        /*Get element on base*/
     2726        Element* basalelement = element->GetBasalElement()->SpawnBasalElement(true);
     2727
     2728        /*Fetch number of nodes and dof for this finite element*/
     2729        int numnodes = basalelement->GetNumberOfNodes();
     2730        int numdof   = numnodes*2;
     2731
     2732        /*Initialize Element matrix and vectors*/
     2733        ElementMatrix* Ke     = basalelement->NewElementMatrix(MLHOApproximationEnum);
     2734        IssmDouble*    B      = xNew<IssmDouble>(3*numdof);
     2735        IssmDouble*    Bprime = xNew<IssmDouble>(3*numdof);
     2736        IssmDouble*    D      = xNewZeroInit<IssmDouble>(3*3);
     2737
     2738        /*Retrieve all inputs and parameters*/
     2739        element->GetVerticesCoordinates(&xyz_list);
     2740        Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input);
     2741        Input* vx_input      = element->GetInput(VxEnum);      _assert_(vx_input);
     2742        Input* vy_input      = element->GetInput(VyEnum);      _assert_(vy_input);
     2743
     2744        /* Start  looping on the number of gaussian points: */
     2745        Gauss* gauss      = element->NewGauss(5);
     2746        Gauss* gauss_base = basalelement->NewGauss();
     2747        while(gauss->next()){
     2748                gauss->SynchronizeGaussBase(gauss_base);
     2749
     2750                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     2751                this->GetBSSA(B,basalelement,2,xyz_list,gauss_base);
     2752                this->GetBSSAprime(Bprime,basalelement,2,xyz_list,gauss_base);
     2753
     2754                element->material->ViscosityL1L2(&viscosity,xyz_list,gauss,vx_input,vy_input,surface_input);
     2755
     2756                for(int i=0;i<3;i++) D[i*3+i]=2*viscosity*gauss->weight*Jdet;
     2757
     2758                TripleMultiply(B,3,numdof,1,
     2759                                        D,3,3,0,
     2760                                        Bprime,3,numdof,0,
     2761                                        &Ke->values[0],1);
     2762        }
     2763
     2764        /*Transform Coordinate System*/
     2765        basalelement->TransformStiffnessMatrixCoord(Ke,XYEnum);
     2766
     2767        /*Clean up and return*/
     2768        delete gauss;
     2769        delete gauss_base;
     2770        basalelement->DeleteMaterials(); delete basalelement;
     2771        xDelete<IssmDouble>(xyz_list);
     2772        xDelete<IssmDouble>(D);
     2773        xDelete<IssmDouble>(Bprime);
     2774        xDelete<IssmDouble>(B);
     2775        return Ke;
     2776}/*}}}*/
     2777ElementVector* StressbalanceAnalysis::CreatePVectorMLHO(Element* element){/*{{{*/
     2778        _error_("not implemented yet");
     2779
     2780        /* Check if ice in element */
     2781        if(!element->IsIceInElement()) return NULL;
     2782
     2783        /*Intermediaries*/
     2784        int      domaintype;
     2785        Element* basalelement;
     2786
     2787        /*Get basal element*/
     2788        element->FindParam(&domaintype,DomainTypeEnum);
     2789        switch(domaintype){
     2790                case Domain2DhorizontalEnum:
     2791                        basalelement = element;
     2792                        break;
     2793                case Domain3DEnum: case Domain2DverticalEnum:
     2794                        if(!element->IsOnBase()) return NULL;
     2795                        basalelement = element->SpawnBasalElement();
     2796                        break;
     2797                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     2798        }
     2799
     2800        /*compute all load vectors for this element*/
     2801        ElementVector* pe1=CreatePVectorMLHODrivingStress(basalelement);
     2802        ElementVector* pe2=CreatePVectorMLHOFront(basalelement);
     2803        ElementVector* pe =new ElementVector(pe1,pe2);
     2804
     2805        /*clean-up and return*/
     2806        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     2807        delete pe1;
     2808        delete pe2;
     2809        return pe;
     2810}/*}}}*/
     2811ElementVector* StressbalanceAnalysis::CreatePVectorMLHODrivingStress(Element* element){/*{{{*/
     2812        _error_("not implemented yet");
     2813
     2814        /*Intermediaries */
     2815        IssmDouble  thickness,Jdet,slope[2];
     2816        IssmDouble* xyz_list = NULL;
     2817
     2818        /*Fetch number of nodes and dof for this finite element*/
     2819        int numnodes = element->GetNumberOfNodes();
     2820
     2821        /*Initialize Element vector and vectors*/
     2822        ElementVector* pe    = element->NewElementVector(MLHOApproximationEnum);
     2823        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     2824
     2825        /*Retrieve all inputs and parameters*/
     2826        element->GetVerticesCoordinates(&xyz_list);
     2827        Input*     thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input);
     2828        Input*     surface_input  =element->GetInput(SurfaceEnum);   _assert_(surface_input);
     2829        IssmDouble rhog = element->FindParam(MaterialsRhoIceEnum)*element->FindParam(ConstantsGEnum);
     2830
     2831        /* Start  looping on the number of gaussian points: */
     2832        Gauss* gauss=element->NewGauss(2);
     2833        while(gauss->next()){
     2834                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     2835                element->NodalFunctions(basis, gauss);
     2836
     2837                thickness_input->GetInputValue(&thickness,gauss);
     2838                surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
     2839
     2840                for(int i=0;i<numnodes;i++){
     2841                        pe->values[i*2+0]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i];
     2842                        pe->values[i*2+1]+=-rhog*thickness*slope[1]*Jdet*gauss->weight*basis[i];
     2843                }
     2844        }
     2845
     2846        /*Transform coordinate system*/
     2847        element->TransformLoadVectorCoord(pe,XYEnum);
     2848
     2849        /*Clean up and return*/
     2850        xDelete<IssmDouble>(xyz_list);
     2851        xDelete<IssmDouble>(basis);
     2852        delete gauss;
     2853        return pe;
     2854}/*}}}*/
     2855ElementVector* StressbalanceAnalysis::CreatePVectorMLHOFront(Element* element){/*{{{*/
     2856        _error_("not implemented yet");
     2857
     2858        /*If no front, return NULL*/
     2859        if(!element->IsIcefront()) return NULL;
     2860
     2861        /*Intermediaries*/
     2862        IssmDouble  Jdet,thickness,bed,sealevel,water_pressure,ice_pressure;
     2863        IssmDouble  surface_under_water,base_under_water,pressure;
     2864        IssmDouble *xyz_list = NULL;
     2865        IssmDouble *xyz_list_front = NULL;
     2866        IssmDouble  normal[2];
     2867
     2868        /*Fetch number of nodes for this finite element*/
     2869        int numnodes = element->GetNumberOfNodes();
     2870
     2871        /*Initialize Element vector and other vectors*/
     2872        ElementVector* pe    = element->NewElementVector(MLHOApproximationEnum);
     2873        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     2874
     2875        /*Retrieve all inputs and parameters*/
     2876        Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
     2877        Input* base_input       = element->GetInput(BaseEnum);       _assert_(base_input);
     2878        Input* sealevel_input       = element->GetInput(SealevelEnum);       _assert_(sealevel_input);
     2879        IssmDouble rho_water   = element->FindParam(MaterialsRhoSeawaterEnum);
     2880        IssmDouble rho_ice     = element->FindParam(MaterialsRhoIceEnum);
     2881        IssmDouble gravity     = element->FindParam(ConstantsGEnum);
     2882        element->GetVerticesCoordinates(&xyz_list);
     2883        element->GetIcefrontCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
     2884        element->NormalSection(&normal[0],xyz_list_front);
     2885
     2886        /*Start looping on Gaussian points*/
     2887        Gauss* gauss=element->NewGauss(xyz_list,xyz_list_front,3);
     2888        while(gauss->next()){
     2889                thickness_input->GetInputValue(&thickness,gauss);
     2890                base_input->GetInputValue(&bed,gauss);
     2891                sealevel_input->GetInputValue(&sealevel,gauss);
     2892                element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
     2893                element->NodalFunctions(basis,gauss);
     2894
     2895                surface_under_water = min(0.,thickness+bed-sealevel); // 0 if the top of the glacier is above water level
     2896                base_under_water    = min(0.,bed-sealevel);           // 0 if the bottom of the glacier is above water level
     2897                water_pressure = 1.0/2.0*gravity*rho_water*(surface_under_water*surface_under_water - base_under_water*base_under_water);
     2898                ice_pressure   = 1.0/2.0*gravity*rho_ice*thickness*thickness;
     2899                pressure = ice_pressure + water_pressure;
     2900
     2901                for (int i=0;i<numnodes;i++){
     2902                        pe->values[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i];
     2903                        pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i];
     2904                }
     2905        }
     2906
     2907        /*Transform coordinate system*/
     2908        element->TransformLoadVectorCoord(pe,XYEnum);
     2909
     2910        /*Clean up and return*/
     2911        xDelete<IssmDouble>(xyz_list);
     2912        xDelete<IssmDouble>(xyz_list_front);
     2913        xDelete<IssmDouble>(basis);
     2914        delete gauss;
     2915        return pe;
     2916}/*}}}*/
     2917void           StressbalanceAnalysis::InputUpdateFromSolutionMLHO(IssmDouble* solution,Element* element){/*{{{*/
     2918
     2919        _error_("not implemented yet");
    25752920
    25762921        int         i,dim,domaintype;

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

    r25514 r25610  
    5656                ElementVector* CreatePVectorL1L2DrivingStress(Element* element);
    5757                void           InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element);
     58                /*MLHO*/
     59                ElementMatrix* CreateKMatrixMLHO(Element* element);
     60                ElementMatrix* CreateKMatrixMLHOFriction(Element* element);
     61                ElementMatrix* CreateKMatrixMLHOViscous(Element* element);
     62                ElementVector* CreatePVectorMLHO(Element* element);
     63                ElementVector* CreatePVectorMLHOFront(Element* element);
     64                ElementVector* CreatePVectorMLHODrivingStress(Element* element);
     65                void           InputUpdateFromSolutionMLHO(IssmDouble* solution,Element* element);
    5866                /*HO*/
    5967                ElementMatrix* CreateJacobianMatrixHO(Element* element);

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

    r25554 r25610  
    1010
    1111        /*Intermediaries*/
    12         bool       isSIA,isSSA,isL1L2,isHO,isFS,iscoupling;
     12        bool       isSIA,isSSA,isL1L2,isMLHO,isHO,isFS,iscoupling;
    1313
    1414        /*Fetch parameters: */
    … …  
    1616        iomodel->FindConstant(&isSSA,"md.flowequation.isSSA");
    1717        iomodel->FindConstant(&isL1L2,"md.flowequation.isL1L2");
     18        iomodel->FindConstant(&isMLHO,"md.flowequation.isMLHO");
    1819        iomodel->FindConstant(&isHO,"md.flowequation.isHO");
    1920        iomodel->FindConstant(&isFS,"md.flowequation.isFS");
    … …  
    2324
    2425        /*Do we have coupling*/
    25         if((isSIA?1.:0.) + (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
     26        if((isSIA?1.:0.) + (isSSA?1.:0.) + (isL1L2?1.:0.) + (isMLHO?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
    2627         iscoupling = true;
    2728        else

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

    r25437 r25610  
    1010
    1111        /*Intermediary*/
    12         bool        isSIA,isSSA,isL1L2,isHO,isFS,iscoupling;
     12        bool        isSIA,isSSA,isL1L2,isMLHO,isHO,isFS,iscoupling;
    1313        int         Mz,Nz;
    1414        IssmDouble *spcvz = NULL;
    … …  
    2121        iomodel->FindConstant(&isSSA,"md.flowequation.isSSA");
    2222        iomodel->FindConstant(&isL1L2,"md.flowequation.isL1L2");
     23        iomodel->FindConstant(&isMLHO,"md.flowequation.isMLHO");
    2324        iomodel->FindConstant(&isHO,"md.flowequation.isHO");
    2425        iomodel->FindConstant(&isFS,"md.flowequation.isFS");
    2526
    2627        /*Do we have coupling*/
    27         if((isSIA?1.:0.) + (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
     28        if((isSIA?1.:0.) + (isSSA?1.:0.) + (isL1L2?1.:0.) + (isMLHO?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
    2829         iscoupling = true;
    2930        else

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r25554 r25610  
    31923192        this->parameters->FindParam(&temp,FlowequationIsSSAEnum); iomodel->AddConstant(new IoConstant(temp,"md.flowequation.isSSA"));
    31933193        this->parameters->FindParam(&temp,FlowequationIsL1L2Enum); iomodel->AddConstant(new IoConstant(temp,"md.flowequation.isL1L2"));
     3194        this->parameters->FindParam(&temp,FlowequationIsMLHOEnum); iomodel->AddConstant(new IoConstant(temp,"md.flowequation.isMLHO"));
    31943195        this->parameters->FindParam(&temp,FlowequationIsHOEnum); iomodel->AddConstant(new IoConstant(temp,"md.flowequation.isHO"));
    31953196        this->parameters->FindParam(&temp,FlowequationIsFSEnum); iomodel->AddConstant(new IoConstant(temp,"md.flowequation.isFS"));

  • issm/trunk-jpl/src/c/classes/Node.cpp

    r25508 r25610  
    117117                        }
    118118                        if(in_approximation==L1L2ApproximationEnum && !reCast<int>(iomodel->Data("md.mesh.vertexonbase")[io_index])){
     119                                this->HardDeactivate();
     120                        }
     121                        if(in_approximation==MLHOApproximationEnum && !reCast<int>(iomodel->Data("md.mesh.vertexonbase")[io_index])){
    119122                                this->HardDeactivate();
    120123                        }
Note: See TracChangeset for help on using the changeset viewer.