Changeset 18227


Ignore:
Timestamp:
07/08/14 07:46:16 (11 years ago)
Author:
seroussi
Message:

NEW: working on Augmented lagrangien sigma nn

Location:
issm/trunk-jpl/src
Files:
2 added
8 edited

Legend:

Unmodified
Added
Removed

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

    r18218 r18227  
    112112        if(fe_FS==XTaylorHoodEnum || fe_FS==LATaylorHoodEnum){
    113113                parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianREnum));
     114                parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianRlambdaEnum));
    114115                parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianThetaEnum));
    115116        }
    … …  
    28792880        /*Intermediaries*/
    28802881        int         i,dim,epssize;
    2881         IssmDouble  r,Jdet,viscosity,DU;
     2882        IssmDouble  r,rl,Jdet,viscosity,DU,DUl;
    28822883        IssmDouble *xyz_list = NULL;
     2884        IssmDouble *xyz_list_base = NULL;
    28832885
    28842886        /*Get problem dimension*/
    … …  
    28932895        if(dim==2) pnumnodes=3;
    28942896        else pnumnodes=6;
     2897        int lnumnodes;
     2898        if(dim==2) lnumnodes=2;
     2899        else lnumnodes=3;
    28952900        //int pnumnodes = element->NumberofNodes(P1Enum);
    28962901        int numdof    = vnumnodes*dim;
    28972902        int pnumdof   = pnumnodes;
     2903        int lnumdof   = lnumnodes;
    28982904
    28992905        /*Prepare coordinate system list*/
    … …  
    29102916        IssmDouble*    BprimeU  = xNew<IssmDouble>(numdof);
    29112917        IssmDouble*    D        = xNewZeroInit<IssmDouble>(epssize*epssize);
     2918        IssmDouble*    CtCUzawa = xNewZeroInit<IssmDouble>(numdof*pnumdof);
     2919        IssmDouble*    C        = xNew<IssmDouble>(lnumdof);
     2920        IssmDouble*    Cprime   = xNew<IssmDouble>(numdof);
    29122921
    29132922        /*Retrieve all inputs and parameters*/
    … …  
    29472956        }
    29482957
     2958        if(element->IsOnBase()){
     2959                element->FindParam(&rl,AugmentedLagrangianRlambdaEnum);
     2960                element->GetVerticesCoordinatesBase(&xyz_list_base);
     2961
     2962                delete gauss;
     2963                Gauss* gauss=element->NewGaussBase(5);
     2964                for(int ig=gauss->begin();ig<gauss->end();ig++){
     2965                        gauss->GaussPoint(ig);
     2966
     2967                        element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
     2968                        this->GetCFS(C,element,dim,xyz_list_base,gauss);
     2969                        this->GetCFSprime(Cprime,element,dim,xyz_list_base,gauss);
     2970                        DUl = gauss->weight*Jdet*sqrt(rl);
     2971                        TripleMultiply(C,1,lnumdof,1,
     2972                                                &DU,1,1,0,
     2973                                                Cprime,1,numdof,0,
     2974                                                CtCUzawa,1);
     2975                }
     2976        }
     2977
    29492978        /*Transform Coordinate System*/
    29502979        element->TransformStiffnessMatrixCoord(Ke,cs_list);
    … …  
    29532982        MatrixMultiply(BtBUzawa,pnumdof,numdof,1,
    29542983                                BtBUzawa,pnumdof,numdof,0,
     2984                                &Ke->values[0],1);
     2985
     2986        /*The sigma naugmentation should not be transformed*/
     2987        MatrixMultiply(CtCUzawa,lnumdof,numdof,1,
     2988                                CtCUzawa,lnumdof,numdof,0,
    29552989                                &Ke->values[0],1);
    29562990
    … …  
    29642998        xDelete<IssmDouble>(BU);
    29652999        xDelete<IssmDouble>(BtBUzawa);
     3000        xDelete<IssmDouble>(Cprime);
     3001        xDelete<IssmDouble>(C);
     3002        xDelete<IssmDouble>(CtCUzawa);
    29663003        xDelete<int>(cs_list);
    29673004        return Ke;
    … …  
    43724409        xDelete<IssmDouble>(vbasis);
    43734410}/*}}}*/
     4411void StressbalanceAnalysis::GetCFS(IssmDouble* C,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     4412        /*Compute C  matrix. C=[Cp1 Cp2 ...] where Cpi=phi_pi.
     4413         */
     4414
     4415        /*Fetch number of nodes for this finite element*/
     4416        int lnumnodes;
     4417        if(dim==2) lnumnodes=3;
     4418        else lnumnodes=6;
     4419        //int pnumnodes = element->NumberofNodes(P1Enum);
     4420
     4421        /*Get nodal functions derivatives*/
     4422        IssmDouble* basis =xNew<IssmDouble>(lnumnodes);
     4423        element->NodalFunctionsP1(basis,gauss);
     4424
     4425        /*Build B: */
     4426        for(int i=0;i<lnumnodes;i++){
     4427                C[i] = basis[i];
     4428        }
     4429
     4430        /*Clean up*/
     4431        xDelete<IssmDouble>(basis);
     4432}/*}}}*/
     4433void StressbalanceAnalysis::GetCFSprime(IssmDouble* Cprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     4434        /*      Compute C'  matrix. C'=[C1' C2' C3']
     4435         *                      Ci' = [  dphi/dx   dphi/dy ]
     4436         *
     4437         *      In 3d
     4438         *         Ci=[ dh/dx   dh/dy    dh/dz  ]
     4439         *      where phi is the finiteelement function for node i.
     4440         */
     4441
     4442        /*Fetch number of nodes for this finite element*/
     4443        int lnumnodes = element->GetNumberOfNodes();
     4444
     4445        /*Get nodal functions derivatives*/
     4446        IssmDouble* dbasis=xNew<IssmDouble>(dim*lnumnodes);
     4447        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     4448
     4449        /*Build C_prime: */
     4450        if(dim==2){
     4451                for(int i=0;i<lnumnodes;i++){
     4452                        Cprime[dim*i+0] = dbasis[0*lnumnodes+i];
     4453                        Cprime[dim*i+1] = dbasis[1*lnumnodes+i];
     4454                }
     4455        }
     4456        else{
     4457                for(int i=0;i<lnumnodes;i++){
     4458                        Cprime[dim*i+0] = dbasis[0*lnumnodes+i];
     4459                        Cprime[dim*i+1] = dbasis[1*lnumnodes+i];
     4460                        Cprime[dim*i+2] = dbasis[2*lnumnodes+i];
     4461                }
     4462        }
     4463
     4464        /*Clean up*/
     4465        xDelete<IssmDouble>(dbasis);
     4466}/*}}}*/
    43744467void StressbalanceAnalysis::GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    43754468

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

    r18184 r18227  
    8686                void GetBFSUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    8787                void GetBFSprimeUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     88                void GetCFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     89                void GetCFSprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    8890                void GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element);
    8991                void InputUpdateFromSolutionFS(IssmDouble* solution,Element* element);

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r18208 r18227  
    669669        AugmentedLagrangianREnum,
    670670        AugmentedLagrangianRhopEnum,
     671        AugmentedLagrangianRlambdaEnum,
     672        AugmentedLagrangianRholambdaEnum,
    671673        AugmentedLagrangianThetaEnum,
    672674        /*}}}*/

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r18208 r18227  
    650650                case AugmentedLagrangianREnum : return "AugmentedLagrangianR";
    651651                case AugmentedLagrangianRhopEnum : return "AugmentedLagrangianRhop";
     652                case AugmentedLagrangianRlambdaEnum : return "AugmentedLagrangianRlambda";
     653                case AugmentedLagrangianRholambdaEnum : return "AugmentedLagrangianRholambda";
    652654                case AugmentedLagrangianThetaEnum : return "AugmentedLagrangianTheta";
    653655                case NoneEnum : return "None";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r18208 r18227  
    665665              else if (strcmp(name,"AugmentedLagrangianR")==0) return AugmentedLagrangianREnum;
    666666              else if (strcmp(name,"AugmentedLagrangianRhop")==0) return AugmentedLagrangianRhopEnum;
     667              else if (strcmp(name,"AugmentedLagrangianRlambda")==0) return AugmentedLagrangianRlambdaEnum;
     668              else if (strcmp(name,"AugmentedLagrangianRholambda")==0) return AugmentedLagrangianRholambdaEnum;
    667669              else if (strcmp(name,"AugmentedLagrangianTheta")==0) return AugmentedLagrangianThetaEnum;
    668670              else if (strcmp(name,"None")==0) return NoneEnum;

  • issm/trunk-jpl/src/m/classes/flowequation.m

    r18208 r18227  
    66classdef flowequation
    77        properties (SetAccess=public)
    8                 isSIA                     = 0;
    9                 isSSA                     = 0;
    10                 isL1L2                    = 0;
    11                 isHO                      = 0;
    12                 isFS                      = 0;
    13                 fe_SSA                    = '';
    14                 fe_HO                     = '';
    15                 fe_FS                     = '';
    16                 augmented_lagrangian_r    = 1.;
    17                 augmented_lagrangian_rhop = 1.;
    18                 XTH_theta                 = 0.;
    19                 vertex_equation           = NaN;
    20                 element_equation          = NaN;
    21                 borderSSA                 = NaN;
    22                 borderHO                  = NaN;
    23                 borderFS                  = NaN;
     8                isSIA                          = 0;
     9                isSSA                          = 0;
     10                isL1L2                         = 0;
     11                isHO                           = 0;
     12                isFS                           = 0;
     13                fe_SSA                         = '';
     14                fe_HO                          = '';
     15                fe_FS                          = '';
     16                augmented_lagrangian_r         = 1.;
     17                augmented_lagrangian_rhop      = 1.;
     18                augmented_lagrangian_rlambda   = 1.;
     19                augmented_lagrangian_rholambda = 1.;
     20                XTH_theta                      = 0.;
     21                vertex_equation                = NaN;
     22                element_equation               = NaN;
     23                borderSSA                      = NaN;
     24                borderHO                       = NaN;
     25                borderFS                       = NaN;
    2426        end
    2527        methods (Static)
    … …  
    137139                                md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',{'P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','LATaylorHood','XTaylorHood','OneLayerP4z','CrouzeixRaviart'});
    138140                                md = checkfield(md,'fieldname','flowequation.augmented_lagrangian_r','numel',[1],'>',0.);
     141                                md = checkfield(md,'fieldname','flowequation.augmented_lagrangian_rlambda','numel',[1],'>',0.);
    139142                                md = checkfield(md,'fieldname','flowequation.augmented_lagrangian_rhop','numel',[1],'>',0.);
     143                                md = checkfield(md,'fieldname','flowequation.augmented_lagrangian_rholambda','numel',[1],'>',0.);
    140144                                md = checkfield(md,'fieldname','flowequation.XTH_theta','numel',[1],'>=',0.,'<',0.5);
    141145                                md = checkfield(md,'fieldname','flowequation.borderSSA','size',[md.mesh.numberofvertices 1],'values',[0 1]);
    … …  
    196200                        WriteData(fid,'enum',AugmentedLagrangianREnum(),'data',obj.augmented_lagrangian_r ,'format','Double');
    197201                        WriteData(fid,'enum',AugmentedLagrangianRhopEnum(),'data',obj.augmented_lagrangian_rhop ,'format','Double');
     202                        WriteData(fid,'enum',AugmentedLagrangianRlambdaEnum(),'data',obj.augmented_lagrangian_rlambda ,'format','Double');
     203                        WriteData(fid,'enum',AugmentedLagrangianRholambdaEnum(),'data',obj.augmented_lagrangian_rholambda ,'format','Double');
    198204                        WriteData(fid,'enum',AugmentedLagrangianThetaEnum() ,'data',obj.XTH_theta ,'format','Double');
    199205                        WriteData(fid,'object',obj,'fieldname','borderSSA','format','DoubleMat','mattype',1);

  • issm/trunk-jpl/src/m/classes/flowequation.py

    r18208 r18227  
    1717        def __init__(self): # {{{
    1818               
    19                 self.isSIA                     = 0
    20                 self.isSSA                     = 0
    21                 self.isL1L2                    = 0
    22                 self.isHO                      = 0
    23                 self.isFS                      = 0
    24                 self.fe_SSA                    = ''
    25                 self.fe_HO                     = ''
    26                 self.fe_FS                     = ''
    27                 self.augmented_lagrangian_r    = 1.
    28                 self.augmented_lagrangian_rhop = 1.
    29                 self.XTH_theta                 = 0.
    30                 self.vertex_equation           = float('NaN')
    31                 self.element_equation          = float('NaN')
    32                 self.borderSSA                 = float('NaN')
    33                 self.borderHO                  = float('NaN')
    34                 self.borderFS                  = float('NaN')
     19                self.isSIA                          = 0
     20                self.isSSA                          = 0
     21                self.isL1L2                         = 0
     22                self.isHO                           = 0
     23                self.isFS                           = 0
     24                self.fe_SSA                         = ''
     25                self.fe_HO                          = ''
     26                self.fe_FS                          = ''
     27                self.augmented_lagrangian_r         = 1.
     28                self.augmented_lagrangian_rhop      = 1.
     29                self.augmented_lagrangian_rlambda   = 1.
     30                self.augmented_lagrangian_rholambda = 1.
     31                self.XTH_theta                      = 0.
     32                self.vertex_equation                = float('NaN')
     33                self.element_equation               = float('NaN')
     34                self.borderSSA                      = float('NaN')
     35                self.borderHO                       = float('NaN')
     36                self.borderFS                       = float('NaN')
    3537
    3638                #set defaults
    … …  
    8587                        md = checkfield(md,'fieldname','flowequation.augmented_lagrangian_r','numel',[1],'>',0.)
    8688                        md = checkfield(md,'fieldname','flowequation.augmented_lagrangian_rhop','numel',[1],'>',0.)
     89                        md = checkfield(md,'fieldname','flowequation.augmented_lagrangian_rlambda','numel',[1],'>',0.)
     90                        md = checkfield(md,'fieldname','flowequation.augmented_lagrangian_rholambda','numel',[1],'>',0.)
    8791                        md = checkfield(md,'fieldname','flowequation.XTH_theta','numel',[1],'>=',0.,'<',.5)
    8892                        if m.strcmp(md.mesh.domaintype(),'2Dhorizontal'):
    … …  
    115119                WriteData(fid,'enum',AugmentedLagrangianREnum(),'data',self.augmented_lagrangian_r ,'format','Double')
    116120                WriteData(fid,'enum',AugmentedLagrangianRhopEnum(),'data',self.augmented_lagrangian_rhop ,'format','Double')
     121                WriteData(fid,'enum',AugmentedLagrangianRlambdaEnum(),'data',self.augmented_lagrangian_rlambda ,'format','Double')
     122                WriteData(fid,'enum',AugmentedLagrangianRholambdaEnum(),'data',self.augmented_lagrangian_rholambda ,'format','Double')
    117123                WriteData(fid,'enum',AugmentedLagrangianThetaEnum() ,'data',self.XTH_theta ,'format','Double')
    118124                WriteData(fid,'object',self,'fieldname','borderSSA','format','DoubleMat','mattype',1)

  • issm/trunk-jpl/src/m/enum/EnumDefinitions.py

    r18208 r18227  
    642642def AugmentedLagrangianREnum(): return StringToEnum("AugmentedLagrangianR")[0]
    643643def AugmentedLagrangianRhopEnum(): return StringToEnum("AugmentedLagrangianRhop")[0]
     644def AugmentedLagrangianRlambdaEnum(): return StringToEnum("AugmentedLagrangianRlambda")[0]
     645def AugmentedLagrangianRholambdaEnum(): return StringToEnum("AugmentedLagrangianRholambda")[0]
    644646def AugmentedLagrangianThetaEnum(): return StringToEnum("AugmentedLagrangianTheta")[0]
    645647def NoneEnum(): return StringToEnum("None")[0]
Note: See TracChangeset for help on using the changeset viewer.