Changeset 17546


Ignore:
Timestamp:
03/26/14 11:35:16 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: done with NewtonSolve

Location:
issm/trunk-jpl/src/c
Files:
1 added
3 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/Makefile.am

    r17525 r17546  
    206206                                        ./shared/Numerics/isnan.cpp\
    207207                                        ./shared/Numerics/cubic.cpp\
     208                                        ./shared/Numerics/NewtonSolveDnorm.cpp\
    208209                                        ./shared/Numerics/extrema.cpp\
    209210                                        ./shared/Numerics/XZvectorsToCoordinateSystem.cpp\

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

    r17541 r17546  
    34463446                sigmapyz_input=element->GetInput(DeviatoricStressyzEnum); _assert_(sigmapyz_input);
    34473447        }
     3448
    34483449        gauss = element->NewGauss();
    34493450        for(int i=0;i<tnumnodes;i++){
    … …  
    42514252
    42524253        /*Intermediaries*/
    4253         int  dim,meshtype;
     4254        int         dim,tausize,meshtype;
     4255        IssmDouble  epsxx,epsyy,epszz,epsxy,epsxz,epsyz,D_scalar;
     4256        IssmDouble  sigmapxx,sigmapyy,sigmapzz,sigmapxy,sigmapxz,sigmapyz;
     4257        IssmDouble  dvx[3],dvy[3],dvz[3],B,n;
     4258        IssmDouble *xyz_list = NULL;
     4259        IssmDouble  Jdet,r;
    42544260
    42554261        parameters->FindParam(&meshtype,MeshTypeEnum);
    42564262        switch(meshtype){
    4257                 case Mesh2DverticalEnum: dim = 2; break;
    4258                 case Mesh3DEnum:         dim = 3; break;
    4259                 case Mesh3DtetrasEnum:   dim = 3; break;
     4263                case Mesh2DverticalEnum: dim = 2; tausize = 3; break;
     4264                case Mesh3DEnum:         dim = 3; tausize = 6; break;
     4265                case Mesh3DtetrasEnum:   dim = 3; tausize = 6; break;
    42604266                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
    42614267        }
     4268
     4269        /*FIXME*/
     4270        r = 1.;
    42624271
    42634272        for(int i=0;i<elements->Size();i++){
    … …  
    42654274
    42664275                /*Get inputs and parameters*/
    4267                 Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input);
    4268                 Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input);
     4276                element->GetVerticesCoordinates(&xyz_list);
     4277                Input*  B_input=element->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
     4278                Input*  n_input=element->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
     4279                Input* vx_input=element->GetInput(VxEnum);                 _assert_(vx_input);
     4280                Input* vy_input=element->GetInput(VyEnum);                 _assert_(vy_input);
    42694281                Input* vz_input;
    42704282                if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
    42714283
    4272                 _error_("NOT implemented yet");
     4284                /*Fetch number of nodes and dof for this finite element*/
     4285                int tnumnodes = element->GetNumberOfVertices();      //Tensors, P1 DG
     4286
     4287                /*Initialize vectors*/
     4288                IssmDouble* tbasis = xNew<IssmDouble>(tnumnodes);
     4289                IssmDouble* Ke     = xNewZeroInit<IssmDouble>(tnumnodes*tnumnodes);
     4290                IssmDouble* pe_xx  = xNewZeroInit<IssmDouble>(tnumnodes);
     4291                IssmDouble* pe_yy  = xNewZeroInit<IssmDouble>(tnumnodes);
     4292                IssmDouble* pe_xy  = xNewZeroInit<IssmDouble>(tnumnodes);
     4293                IssmDouble* pe_zz  = NULL; IssmDouble* pe_xz  = NULL; IssmDouble* pe_yz  = NULL;
     4294                if(dim==3){
     4295                        pe_zz = xNewZeroInit<IssmDouble>(tnumnodes);
     4296                        pe_xz = xNewZeroInit<IssmDouble>(tnumnodes);
     4297                        pe_yz = xNewZeroInit<IssmDouble>(tnumnodes);
     4298                }
     4299
     4300                /*Get tau*/
     4301                Input* sigmapxx_input=element->GetInput(DeviatoricStressxxEnum); _assert_(sigmapxx_input);
     4302                Input* sigmapyy_input=element->GetInput(DeviatoricStressyyEnum); _assert_(sigmapyy_input);
     4303                Input* sigmapxy_input=element->GetInput(DeviatoricStressxyEnum); _assert_(sigmapxy_input);
     4304                Input* sigmapzz_input=NULL; Input* sigmapxz_input=NULL; Input* sigmapyz_input=NULL;
     4305                if(dim==3){
     4306                        sigmapzz_input=element->GetInput(DeviatoricStresszzEnum); _assert_(sigmapzz_input);
     4307                        sigmapxz_input=element->GetInput(DeviatoricStressxzEnum); _assert_(sigmapxz_input);
     4308                        sigmapyz_input=element->GetInput(DeviatoricStressyzEnum); _assert_(sigmapyz_input);
     4309                }
     4310
     4311                Gauss* gauss=element->NewGauss(5);
     4312                for(int ig=gauss->begin();ig<gauss->end();ig++){
     4313                        gauss->GaussPoint(ig);
     4314                        element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     4315                        element->NodalFunctionsTensor(tbasis,gauss);
     4316
     4317                        /*Get tau from inputs*/
     4318                        sigmapxx_input->GetInputValue(&sigmapxx,gauss);
     4319                        sigmapyy_input->GetInputValue(&sigmapyy,gauss);
     4320                        sigmapxy_input->GetInputValue(&sigmapxy,gauss);
     4321                        if(dim==3){
     4322                                sigmapzz_input->GetInputValue(&sigmapzz,gauss);
     4323                                sigmapxz_input->GetInputValue(&sigmapxz,gauss);
     4324                                sigmapyz_input->GetInputValue(&sigmapyz,gauss);
     4325                        }
     4326
     4327                        /*Calculate d from previous results*/
     4328                        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
     4329                        vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
     4330                        if(dim==3){
     4331                                vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
     4332                        }
     4333                        epsxx = dvx[0];
     4334                        epsyy = dvy[1];
     4335                        epsxy = 0.5*(dvx[1] + dvy[0]);
     4336                        if(dim==3){
     4337                                epszz = dvz[2];               
     4338                                epsxz = 0.5*(dvx[2] + dvz[0]);
     4339                                epsyz = 0.5*(dvy[2] + dvz[1]);
     4340                        }
     4341
     4342                        /*Solve 2 eta_0 |d|^s-1 + r |d| = |rD(u) + tau|*/
     4343                        IssmDouble coef1,coef2,coef3;
     4344                        B_input->GetInputValue(&B,gauss);
     4345                        n_input->GetInputValue(&n,gauss);
     4346                        coef1 = (B/2.)*pow(1./sqrt(2.),(1.-n)/n); //2 eta_0
     4347                        coef2 = r;
     4348                        if(dim==2){
     4349                                coef3 = sqrt(
     4350                                                          (r*epsxx + sigmapxx)*(r*epsxx + sigmapxx)
     4351                                                        + (r*epsyy + sigmapyy)*(r*epsyy + sigmapyy)
     4352                                                        + 2*(r*epsxy + sigmapxy)*(r*epsxy + sigmapxy)
     4353                                                        );
     4354                        }
     4355                        else{
     4356                                coef3 = sqrt(
     4357                                                          (r*epsxx + sigmapxx)*(r*epsxx + sigmapxx)
     4358                                                        + (r*epsyy + sigmapyy)*(r*epsyy + sigmapyy)
     4359                                                        + (r*epszz + sigmapzz)*(r*epszz + sigmapzz)
     4360                                                        + 2*(r*epsxy + sigmapxy)*(r*epsxy + sigmapxy)
     4361                                                        + 2*(r*epsxz + sigmapxz)*(r*epsxz + sigmapxz)
     4362                                                        + 2*(r*epsyz + sigmapyz)*(r*epsyz + sigmapyz)
     4363                                                        );
     4364                        }
     4365                        IssmDouble dnorm;
     4366                        NewtonSolveDnorm(&dnorm,coef1,coef2,coef3,n);
     4367
     4368                        /*Create Ke*/
     4369                        D_scalar=(coef1*pow(dnorm,(1.-n)/n)+r)*gauss->weight*Jdet;
     4370                        TripleMultiply(tbasis,tnumnodes,1,0,
     4371                                                &D_scalar,1,1,0,
     4372                                                tbasis,1,tnumnodes,0,
     4373                                                Ke,1);
     4374
     4375                        /*Create Right hand sides*/
     4376                        for(int ii=0;ii<tnumnodes;ii++) pe_xx[ii] += sigmapxx*tbasis[ii]*gauss->weight*Jdet;
     4377                        for(int ii=0;ii<tnumnodes;ii++) pe_yy[ii] += sigmapyy*tbasis[ii]*gauss->weight*Jdet;
     4378                        for(int ii=0;ii<tnumnodes;ii++) pe_xy[ii] += sigmapxy*tbasis[ii]*gauss->weight*Jdet;
     4379                        if(dim==3){
     4380                                for(int ii=0;ii<tnumnodes;ii++) pe_zz[ii] += sigmapzz*tbasis[ii]*gauss->weight*Jdet;
     4381                                for(int ii=0;ii<tnumnodes;ii++) pe_xz[ii] += sigmapxz*tbasis[ii]*gauss->weight*Jdet;
     4382                                for(int ii=0;ii<tnumnodes;ii++) pe_yz[ii] += sigmapyz*tbasis[ii]*gauss->weight*Jdet;
     4383                        }
     4384                }
     4385
     4386                /*Solve the systems*/
     4387                _error_("S");
     4388
     4389                /*Clean up and */
     4390                xDelete<IssmDouble>(tbasis);
     4391                xDelete<IssmDouble>(Ke);
     4392                xDelete<IssmDouble>(pe_xx);
     4393                xDelete<IssmDouble>(pe_yy);
     4394                xDelete<IssmDouble>(pe_zz);
     4395                xDelete<IssmDouble>(pe_xy);
     4396                xDelete<IssmDouble>(pe_xz);
     4397                xDelete<IssmDouble>(pe_yz);
    42734398        }
    42744399

  • issm/trunk-jpl/src/c/shared/Numerics/numerics.h

    r15128 r17546  
    3838int         cubic(IssmDouble a, IssmDouble b, IssmDouble c, IssmDouble d, double X[3], int *num);
    3939
     40int         NewtonSolveDnorm(IssmDouble* pdnorm,IssmDouble c1,IssmDouble c2,IssmDouble c3,IssmDouble n);
     41
    4042#endif //ifndef _NUMERICS_H_
Note: See TracChangeset for help on using the changeset viewer.