Changeset 26989


Ignore:
Timestamp:
05/05/22 05:52:29 (3 years ago)
Author:
Mathieu Morlighem
Message:

BUG: fixed memory leaks

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r26478 r26989  
    226226
    227227        /*Intermediaries */
    228         bool        incomplete_adjoint;
    229228        IssmDouble  Jdet,mu_prime,n,thickness,mu,effmu;
    230229        IssmDouble *xyz_list = NULL;
    … …  
    238237        IssmDouble  vxshear, vyshear;
    239238
    240    Element* basalelement;
     239   /*Get basal element*/
     240        Element* basalelement = NULL;
     241   element->FindParam(&domaintype,DomainTypeEnum);
     242   switch(domaintype){
     243      case Domain2DhorizontalEnum:
     244         basalelement = element;
     245         break;
     246      case Domain3DEnum: case Domain2DverticalEnum:
     247         _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     248         break;
     249      default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     250   }
     251
     252        /*Initialize Jacobian with regular HO (first part of the Gateau derivative)*/
     253        bool        incomplete_adjoint;
     254        element->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
     255        StressbalanceAnalysis* analysis = new StressbalanceAnalysis();
     256        ElementMatrix* Ke=analysis->CreateKMatrix(basalelement);
     257        delete analysis;
     258        if(incomplete_adjoint) return Ke;
     259
     260        /*Second part of the Gateau derivative if requested*/
     261
     262        /*Fetch number of nodes and dof for this finite element*/
     263   int numnodes = basalelement->GetNumberOfNodes();
     264   IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes); // like SSA
     265   IssmDouble* basis  = xNew<IssmDouble>(numnodes); // like SSA
     266
     267        /*Retrieve all inputs and parameters*/
     268        element->GetVerticesCoordinates(&xyz_list);
     269   Input* vx_input       = element->GetInput(VxEnum);        _assert_(vx_input); //vertically integrated vx
     270   Input* vy_input       = element->GetInput(VyEnum);        _assert_(vy_input); //vertically integrated vy
     271   Input* vxbase_input   = element->GetInput(VxBaseEnum);    _assert_(vxbase_input);
     272   Input* vybase_input   = element->GetInput(VyBaseEnum);    _assert_(vybase_input);
     273   Input* vxshear_input  = element->GetInput(VxShearEnum);   _assert_(vxshear_input);
     274   Input* vyshear_input  = element->GetInput(VyShearEnum);   _assert_(vyshear_input);
     275   Input* thickness_input= element->GetInput(ThicknessEnum); _assert_(thickness_input);
     276   Input* n_input        = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
     277
     278        /* Start  looping on the number of gaussian points: */
     279   Gauss* gauss      = element->NewGauss(5);
     280   Gauss* gauss_base = basalelement->NewGauss();
     281   while(gauss->next()){
     282      gauss->SynchronizeGaussBase(gauss_base);
     283
     284      element->JacobianDeterminant(&Jdet,xyz_list,gauss_base);
     285      basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss_base);
     286      basalelement->NodalFunctions(basis, gauss_base);
     287
     288      thickness_input->GetInputValue(&thickness, gauss);
     289      n_input->GetInputValue(&n,gauss);
     290
     291      vxshear_input->GetInputValue(&vxshear,gauss);
     292      vyshear_input->GetInputValue(&vyshear,gauss);
     293
     294                element->material->ViscosityMLHOAdjoint(&viscosity[0],dim,xyz_list,gauss,vxbase_input,vybase_input,vxshear_input,vyshear_input,thickness_input,n_input);
     295
     296                effmu = 2.0*(1-n)/2.0/n;
     297
     298                element->StrainRateHO(&epsilonbase[0],xyz_list,gauss,vxbase_input, vybase_input);
     299                element->StrainRateHO(&epsilonshear[0],xyz_list,gauss,vxshear_input, vyshear_input);
     300
     301                e1b[0] = 2*epsilonbase[0]+epsilonbase[1];               e1b[1] = epsilonbase[2];
     302                e2b[1] = epsilonbase[0]+2*epsilonbase[1];               e2b[0] = epsilonbase[2];
     303
     304                e1sh[0] = 2*epsilonshear[0]+epsilonshear[1];            e1sh[1] = epsilonshear[2];
     305                e2sh[1] = epsilonshear[0]+2*epsilonshear[1];            e2sh[0] = epsilonshear[2];
     306
     307                for(int i=0;i<numnodes;i++){
     308                        for(int j=0;j<numnodes;j++){
     309                                eb1i = e1b[0]*dbasis[0*numnodes+i]+e1b[1]*dbasis[1*numnodes+i];
     310                                eb1j = e1b[0]*dbasis[0*numnodes+j]+e1b[1]*dbasis[1*numnodes+j];
     311                                eb2i = e2b[0]*dbasis[0*numnodes+i]+e2b[1]*dbasis[1*numnodes+i];
     312                                eb2j = e2b[0]*dbasis[0*numnodes+j]+e2b[1]*dbasis[1*numnodes+j];
     313            esh1i = e1sh[0]*dbasis[0*numnodes+i]+e1sh[1]*dbasis[1*numnodes+i];
     314            esh1j = e1sh[0]*dbasis[0*numnodes+j]+e1sh[1]*dbasis[1*numnodes+j];
     315            esh2i = e2sh[0]*dbasis[0*numnodes+i]+e2sh[1]*dbasis[1*numnodes+i];
     316            esh2j = e2sh[0]*dbasis[0*numnodes+j]+e2sh[1]*dbasis[1*numnodes+j];
     317
     318                                Ke->values[4*numnodes*(4*i+0)+4*j+0]+=gauss->weight*Jdet*effmu*(viscosity[0]*eb1j*eb1i+viscosity[1]*(esh1j*eb1i+eb1j*esh1i)+viscosity[2]*esh1j*esh1i);
     319                                Ke->values[4*numnodes*(4*i+1)+4*j+0]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb1j*eb1i+viscosity[2]*(esh1j*eb1i+eb1j*esh1i)+viscosity[4]*esh1j*esh1i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*eb1j*basis[i]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*esh1j*basis[i]);
     320                                Ke->values[4*numnodes*(4*i+2)+4*j+0]+=gauss->weight*Jdet*effmu*(viscosity[0]*eb1j*eb2i+viscosity[1]*(esh1j*eb2i+eb1j*esh2i)+viscosity[2]*esh1j*esh2i);
     321                                Ke->values[4*numnodes*(4*i+3)+4*j+0]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb1j*eb2i+viscosity[2]*(esh1j*eb2i+eb1j*esh2i)+viscosity[4]*esh1j*esh2i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*eb1j*basis[i]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*esh1j*basis[i]);
     322
     323                                Ke->values[4*numnodes*(4*i+0)+4*j+1]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb1j*eb1i+viscosity[2]*(esh1j*eb1i+eb1j*esh1i)+viscosity[4]*esh1j*esh1i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*eb1i*basis[j]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*esh1i*basis[j]);
     324                                Ke->values[4*numnodes*(4*i+1)+4*j+1]+=gauss->weight*Jdet*effmu*(viscosity[2]*eb1j*eb1i+viscosity[4]*(esh1j*eb1i+eb1j*esh1i)+viscosity[5]*esh1j*esh1i+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*(eb1i*basis[j]+eb1j*basis[i])+viscosity[7]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*(esh1i*basis[j]+esh1j*basis[i])+viscosity[8]*(n+1)*(n+1)*(n+1)*(n+1)/4.0/thickness/thickness/thickness/thickness*vxshear*vxshear*basis[j]*basis[i]);
     325                                Ke->values[4*numnodes*(4*i+2)+4*j+1]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb1j*eb2i+viscosity[2]*(esh1j*eb2i+eb1j*esh2i)+viscosity[4]*esh1j*esh2i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*eb2i*basis[j]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*esh2i*basis[j]);
     326                                Ke->values[4*numnodes*(4*i+3)+4*j+1]+=gauss->weight*Jdet*effmu*(viscosity[2]*eb1j*eb2i+viscosity[4]*(esh1j*eb2i+eb1j*esh2i)+viscosity[5]*esh1j*esh2i+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*(vxshear*eb2i*basis[j]+vyshear*eb1j*basis[i])+viscosity[7]*(n+1)*(n+1)/2.0/thickness/thickness*(vxshear*esh2i*basis[j]+vyshear*esh1j*basis[i])+viscosity[8]*(n+1)*(n+1)*(n+1)*(n+1)/4.0/thickness/thickness/thickness/thickness*vxshear*vyshear*basis[j]*basis[i]);
     327                               
     328
     329                                Ke->values[4*numnodes*(4*i+0)+4*j+2]+=gauss->weight*Jdet*effmu*(viscosity[0]*eb2j*eb1i+viscosity[1]*(esh2j*eb1i+eb2j*esh1i)+viscosity[2]*esh2j*esh1i);
     330                                Ke->values[4*numnodes*(4*i+1)+4*j+2]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb2j*eb1i+viscosity[2]*(esh2j*eb1i+eb2j*esh1i)+viscosity[4]*esh2j*esh1i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*eb2j*basis[i]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*esh2j*basis[i]);
     331                                Ke->values[4*numnodes*(4*i+2)+4*j+2]+=gauss->weight*Jdet*effmu*(viscosity[0]*eb2j*eb2i+viscosity[1]*(esh2j*eb2i+eb2j*esh2i)+viscosity[2]*esh2j*esh2i);
     332                                Ke->values[4*numnodes*(4*i+3)+4*j+2]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb2j*eb2i+viscosity[2]*(esh2j*eb2i+eb2j*esh2i)+viscosity[4]*esh2j*esh2i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*eb2j*basis[i]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*esh2j*basis[i]);
     333
     334
     335                                Ke->values[4*numnodes*(4*i+0)+4*j+3]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb2j*eb1i+viscosity[2]*(esh2j*eb1i+eb2j*esh1i)+viscosity[4]*esh2j*esh1i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*eb1i*basis[j]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*esh1i*basis[j]);
     336                                Ke->values[4*numnodes*(4*i+1)+4*j+3]+=gauss->weight*Jdet*effmu*(viscosity[2]*eb2j*eb1i+viscosity[4]*(esh2j*eb1i+eb2j*esh1i)+viscosity[5]*esh2j*esh1i+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*(vyshear*eb1i*basis[j]+vxshear*eb2j*basis[i])+viscosity[7]*(n+1)*(n+1)/2.0/thickness/thickness*(vxshear*esh2j*basis[i]+vyshear*esh1i*basis[j])+viscosity[8]*(n+1)*(n+1)*(n+1)*(n+1)/4.0/thickness/thickness/thickness/thickness*vxshear*vxshear*basis[j]*basis[i]);
     337                                Ke->values[4*numnodes*(4*i+2)+4*j+3]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb2j*eb2i+viscosity[2]*(esh2j*eb2i+eb2j*esh2i)+viscosity[4]*esh2j*esh2i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*eb2i*basis[j]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*esh2i*basis[j]);
     338                                Ke->values[4*numnodes*(4*i+3)+4*j+3]+=gauss->weight*Jdet*effmu*(viscosity[2]*eb2j*eb2i+viscosity[4]*(esh2j*eb2i+eb2j*esh2i)+viscosity[5]*esh2j*esh2i+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*(eb2i*basis[j]+eb1j*basis[i])+viscosity[7]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*(esh2i*basis[j]+esh1j*basis[i])+viscosity[8]*(n+1)*(n+1)*(n+1)*(n+1)/4.0/thickness/thickness/thickness/thickness*vyshear*vyshear*basis[j]*basis[i]);
     339                        }
     340                }
     341        }
     342
     343        /*Transform Coordinate System*/
     344        /*element->TransformStiffnessMatrixCoord(Ke,XYEnum);*/
     345
     346   /*Clean up and return*/
     347   delete gauss;
     348   delete gauss_base;
     349   if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;};
     350   xDelete<IssmDouble>(xyz_list);
     351   xDelete<IssmDouble>(dbasis);
     352   xDelete<IssmDouble>(basis);
     353   return Ke;
     354}/*}}}*/
     355ElementMatrix* AdjointHorizAnalysis::CreateKMatrixMLHOVerticalIntergrated(Element* element){/*{{{*/
     356
     357        /* Check if ice in element */
     358        if(!element->IsIceInElement()) return NULL;
     359
     360        /*Intermediaries */
     361        IssmDouble  Jdet,mu_prime,n,thickness,mu,effmu;
     362        IssmDouble *xyz_list = NULL;
     363        IssmDouble      zeta, epsilon_eff;
     364   int                  domaintype;
     365        int                     dim=2;
     366
     367        IssmDouble  e1phi1i, e1phi1j, e2phi1i, e2phi1j, e1phi2i, e1phi2j, e2phi2i, e2phi2j;
     368        IssmDouble  epsilon[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/
     369        IssmDouble      e1[3],e2[3];
     370
     371   Element* basalelement = NULL;
    241372
    242373   /*Get basal element*/
    … …  
    251382      default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    252383   }
    253         /*Fetch number of nodes and dof for this finite element*/
    254    int numnodes = basalelement->GetNumberOfNodes();
    255    IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes); // like SSA
    256    IssmDouble* basis  = xNew<IssmDouble>(numnodes); // like SSA
    257384
    258385        /*Initialize Jacobian with regular HO (first part of the Gateau derivative)*/
     386        bool incomplete_adjoint;
    259387        element->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
    260388        StressbalanceAnalysis* analysis = new StressbalanceAnalysis();
    … …  
    263391        if(incomplete_adjoint) return Ke;
    264392
    265         /*Retrieve all inputs and parameters*/
    266         element->GetVerticesCoordinates(&xyz_list);
    267    Input* vx_input       = element->GetInput(VxEnum);        _assert_(vx_input); //vertically integrated vx
    268    Input* vy_input       = element->GetInput(VyEnum);        _assert_(vy_input); //vertically integrated vy
    269    Input* vxbase_input   = element->GetInput(VxBaseEnum);    _assert_(vxbase_input);
    270    Input* vybase_input   = element->GetInput(VyBaseEnum);    _assert_(vybase_input);
    271    Input* vxshear_input  = element->GetInput(VxShearEnum);   _assert_(vxshear_input);
    272    Input* vyshear_input  = element->GetInput(VyShearEnum);   _assert_(vyshear_input);
    273    Input* thickness_input= element->GetInput(ThicknessEnum); _assert_(thickness_input);
    274    Input* n_input        = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
    275 
    276         /* Start  looping on the number of gaussian points: */
    277    Gauss* gauss      = element->NewGauss(5);
    278    Gauss* gauss_base = basalelement->NewGauss();
    279    while(gauss->next()){
    280       gauss->SynchronizeGaussBase(gauss_base);
    281 
    282       element->JacobianDeterminant(&Jdet,xyz_list,gauss_base);
    283       basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss_base);
    284       basalelement->NodalFunctions(basis, gauss_base);
    285 
    286       thickness_input->GetInputValue(&thickness, gauss);
    287       n_input->GetInputValue(&n,gauss);
    288 
    289       vxshear_input->GetInputValue(&vxshear,gauss);
    290       vyshear_input->GetInputValue(&vyshear,gauss);
    291 
    292                 element->material->ViscosityMLHOAdjoint(&viscosity[0],dim,xyz_list,gauss,vxbase_input,vybase_input,vxshear_input,vyshear_input,thickness_input,n_input);
    293 
    294                 effmu = 2.0*(1-n)/2.0/n;
    295 
    296                 element->StrainRateHO(&epsilonbase[0],xyz_list,gauss,vxbase_input, vybase_input);
    297                 element->StrainRateHO(&epsilonshear[0],xyz_list,gauss,vxshear_input, vyshear_input);
    298 
    299                 e1b[0] = 2*epsilonbase[0]+epsilonbase[1];               e1b[1] = epsilonbase[2];
    300                 e2b[1] = epsilonbase[0]+2*epsilonbase[1];               e2b[0] = epsilonbase[2];
    301 
    302                 e1sh[0] = 2*epsilonshear[0]+epsilonshear[1];            e1sh[1] = epsilonshear[2];
    303                 e2sh[1] = epsilonshear[0]+2*epsilonshear[1];            e2sh[0] = epsilonshear[2];
    304 
    305                 for(int i=0;i<numnodes;i++){
    306                         for(int j=0;j<numnodes;j++){
    307                                 eb1i = e1b[0]*dbasis[0*numnodes+i]+e1b[1]*dbasis[1*numnodes+i];
    308                                 eb1j = e1b[0]*dbasis[0*numnodes+j]+e1b[1]*dbasis[1*numnodes+j];
    309                                 eb2i = e2b[0]*dbasis[0*numnodes+i]+e2b[1]*dbasis[1*numnodes+i];
    310                                 eb2j = e2b[0]*dbasis[0*numnodes+j]+e2b[1]*dbasis[1*numnodes+j];
    311             esh1i = e1sh[0]*dbasis[0*numnodes+i]+e1sh[1]*dbasis[1*numnodes+i];
    312             esh1j = e1sh[0]*dbasis[0*numnodes+j]+e1sh[1]*dbasis[1*numnodes+j];
    313             esh2i = e2sh[0]*dbasis[0*numnodes+i]+e2sh[1]*dbasis[1*numnodes+i];
    314             esh2j = e2sh[0]*dbasis[0*numnodes+j]+e2sh[1]*dbasis[1*numnodes+j];
    315 
    316                                 Ke->values[4*numnodes*(4*i+0)+4*j+0]+=gauss->weight*Jdet*effmu*(viscosity[0]*eb1j*eb1i+viscosity[1]*(esh1j*eb1i+eb1j*esh1i)+viscosity[2]*esh1j*esh1i);
    317                                 Ke->values[4*numnodes*(4*i+1)+4*j+0]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb1j*eb1i+viscosity[2]*(esh1j*eb1i+eb1j*esh1i)+viscosity[4]*esh1j*esh1i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*eb1j*basis[i]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*esh1j*basis[i]);
    318                                 Ke->values[4*numnodes*(4*i+2)+4*j+0]+=gauss->weight*Jdet*effmu*(viscosity[0]*eb1j*eb2i+viscosity[1]*(esh1j*eb2i+eb1j*esh2i)+viscosity[2]*esh1j*esh2i);
    319                                 Ke->values[4*numnodes*(4*i+3)+4*j+0]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb1j*eb2i+viscosity[2]*(esh1j*eb2i+eb1j*esh2i)+viscosity[4]*esh1j*esh2i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*eb1j*basis[i]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*esh1j*basis[i]);
    320 
    321                                 Ke->values[4*numnodes*(4*i+0)+4*j+1]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb1j*eb1i+viscosity[2]*(esh1j*eb1i+eb1j*esh1i)+viscosity[4]*esh1j*esh1i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*eb1i*basis[j]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*esh1i*basis[j]);
    322                                 Ke->values[4*numnodes*(4*i+1)+4*j+1]+=gauss->weight*Jdet*effmu*(viscosity[2]*eb1j*eb1i+viscosity[4]*(esh1j*eb1i+eb1j*esh1i)+viscosity[5]*esh1j*esh1i+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*(eb1i*basis[j]+eb1j*basis[i])+viscosity[7]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*(esh1i*basis[j]+esh1j*basis[i])+viscosity[8]*(n+1)*(n+1)*(n+1)*(n+1)/4.0/thickness/thickness/thickness/thickness*vxshear*vxshear*basis[j]*basis[i]);
    323                                 Ke->values[4*numnodes*(4*i+2)+4*j+1]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb1j*eb2i+viscosity[2]*(esh1j*eb2i+eb1j*esh2i)+viscosity[4]*esh1j*esh2i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*eb2i*basis[j]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*esh2i*basis[j]);
    324                                 Ke->values[4*numnodes*(4*i+3)+4*j+1]+=gauss->weight*Jdet*effmu*(viscosity[2]*eb1j*eb2i+viscosity[4]*(esh1j*eb2i+eb1j*esh2i)+viscosity[5]*esh1j*esh2i+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*(vxshear*eb2i*basis[j]+vyshear*eb1j*basis[i])+viscosity[7]*(n+1)*(n+1)/2.0/thickness/thickness*(vxshear*esh2i*basis[j]+vyshear*esh1j*basis[i])+viscosity[8]*(n+1)*(n+1)*(n+1)*(n+1)/4.0/thickness/thickness/thickness/thickness*vxshear*vyshear*basis[j]*basis[i]);
    325                                
    326 
    327                                 Ke->values[4*numnodes*(4*i+0)+4*j+2]+=gauss->weight*Jdet*effmu*(viscosity[0]*eb2j*eb1i+viscosity[1]*(esh2j*eb1i+eb2j*esh1i)+viscosity[2]*esh2j*esh1i);
    328                                 Ke->values[4*numnodes*(4*i+1)+4*j+2]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb2j*eb1i+viscosity[2]*(esh2j*eb1i+eb2j*esh1i)+viscosity[4]*esh2j*esh1i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*eb2j*basis[i]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*esh2j*basis[i]);
    329                                 Ke->values[4*numnodes*(4*i+2)+4*j+2]+=gauss->weight*Jdet*effmu*(viscosity[0]*eb2j*eb2i+viscosity[1]*(esh2j*eb2i+eb2j*esh2i)+viscosity[2]*esh2j*esh2i);
    330                                 Ke->values[4*numnodes*(4*i+3)+4*j+2]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb2j*eb2i+viscosity[2]*(esh2j*eb2i+eb2j*esh2i)+viscosity[4]*esh2j*esh2i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*eb2j*basis[i]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*esh2j*basis[i]);
    331 
    332 
    333                                 Ke->values[4*numnodes*(4*i+0)+4*j+3]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb2j*eb1i+viscosity[2]*(esh2j*eb1i+eb2j*esh1i)+viscosity[4]*esh2j*esh1i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*eb1i*basis[j]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*esh1i*basis[j]);
    334                                 Ke->values[4*numnodes*(4*i+1)+4*j+3]+=gauss->weight*Jdet*effmu*(viscosity[2]*eb2j*eb1i+viscosity[4]*(esh2j*eb1i+eb2j*esh1i)+viscosity[5]*esh2j*esh1i+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*(vyshear*eb1i*basis[j]+vxshear*eb2j*basis[i])+viscosity[7]*(n+1)*(n+1)/2.0/thickness/thickness*(vxshear*esh2j*basis[i]+vyshear*esh1i*basis[j])+viscosity[8]*(n+1)*(n+1)*(n+1)*(n+1)/4.0/thickness/thickness/thickness/thickness*vxshear*vxshear*basis[j]*basis[i]);
    335                                 Ke->values[4*numnodes*(4*i+2)+4*j+3]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb2j*eb2i+viscosity[2]*(esh2j*eb2i+eb2j*esh2i)+viscosity[4]*esh2j*esh2i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*eb2i*basis[j]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*esh2i*basis[j]);
    336                                 Ke->values[4*numnodes*(4*i+3)+4*j+3]+=gauss->weight*Jdet*effmu*(viscosity[2]*eb2j*eb2i+viscosity[4]*(esh2j*eb2i+eb2j*esh2i)+viscosity[5]*esh2j*esh2i+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*(eb2i*basis[j]+eb1j*basis[i])+viscosity[7]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*(esh2i*basis[j]+esh1j*basis[i])+viscosity[8]*(n+1)*(n+1)*(n+1)*(n+1)/4.0/thickness/thickness/thickness/thickness*vyshear*vyshear*basis[j]*basis[i]);
    337                         }
    338                 }
    339         }
    340 
    341         /*Transform Coordinate System*/
    342         /*element->TransformStiffnessMatrixCoord(Ke,XYEnum);*/
    343 
    344    /*Clean up and return*/
    345    delete gauss;
    346    delete gauss_base;
    347    if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;};
    348    xDelete<IssmDouble>(xyz_list);
    349    xDelete<IssmDouble>(dbasis);
    350    xDelete<IssmDouble>(basis);
    351    return Ke;
    352 }/*}}}*/
    353 ElementMatrix* AdjointHorizAnalysis::CreateKMatrixMLHOVerticalIntergrated(Element* element){/*{{{*/
    354 
    355         /* Check if ice in element */
    356         if(!element->IsIceInElement()) return NULL;
    357 
    358         /*Intermediaries */
    359         bool        incomplete_adjoint;
    360         IssmDouble  Jdet,mu_prime,n,thickness,mu,effmu;
    361         IssmDouble *xyz_list = NULL;
    362         IssmDouble      zeta, epsilon_eff;
    363    int                  domaintype;
    364         int                     dim=2;
    365 
    366         IssmDouble  e1phi1i, e1phi1j, e2phi1i, e2phi1j, e1phi2i, e1phi2j, e2phi2i, e2phi2j;
    367         IssmDouble  epsilon[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/
    368         IssmDouble      e1[3],e2[3];
    369 
    370    Element* basalelement;
    371 
    372    /*Get basal element*/
    373    element->FindParam(&domaintype,DomainTypeEnum);
    374    switch(domaintype){
    375       case Domain2DhorizontalEnum:
    376          basalelement = element;
    377          break;
    378       case Domain3DEnum: case Domain2DverticalEnum:
    379          _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    380          break;
    381       default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    382    }
     393        /*Second part of the Gateau derivative if requested*/
     394
     395
    383396        /*Fetch number of nodes and dof for this finite element*/
    384397   int numnodes = basalelement->GetNumberOfNodes();
    385398   IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes); // like SSA
    386399   IssmDouble* basis  = xNew<IssmDouble>(numnodes); // like SSA
    387 
    388         /*Initialize Jacobian with regular HO (first part of the Gateau derivative)*/
    389         element->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
    390         StressbalanceAnalysis* analysis = new StressbalanceAnalysis();
    391         ElementMatrix* Ke=analysis->CreateKMatrix(basalelement);
    392         delete analysis;
    393         if(incomplete_adjoint) return Ke;
    394400
    395401        /*Retrieve all inputs and parameters*/
    … …  
    34283434        /*Fetch dof list and allocate solution vectors*/
    34293435        element->GetDofListLocal(&doflist,MLHOApproximationEnum,GsetEnum);
    3430         IssmDouble* values  = xNew<IssmDouble>(numdof);
    3431         IssmDouble* lambdax = xNew<IssmDouble>(numnodes);
    3432         IssmDouble* lambday = xNew<IssmDouble>(numnodes);
    3433         IssmDouble* lambdabx = xNew<IssmDouble>(numnodes);
    3434         IssmDouble* lambdaby = xNew<IssmDouble>(numnodes);
     3436        IssmDouble* values    = xNew<IssmDouble>(numdof);
     3437        IssmDouble* lambdax   = xNew<IssmDouble>(numnodes);
     3438        IssmDouble* lambday   = xNew<IssmDouble>(numnodes);
     3439        IssmDouble* lambdabx  = xNew<IssmDouble>(numnodes);
     3440        IssmDouble* lambdaby  = xNew<IssmDouble>(numnodes);
    34353441        IssmDouble* lambdashx = xNew<IssmDouble>(numnodes);
    34363442        IssmDouble* lambdashy = xNew<IssmDouble>(numnodes);
    3437         IssmDouble* n                    = xNew<IssmDouble>(numnodes);
     3443        IssmDouble* n        = xNew<IssmDouble>(numnodes);
    34383444
    34393445        /*Use the dof list to index into the solution vector: */
    … …  
    34853491        xDelete<IssmDouble>(lambdashx);
    34863492        xDelete<IssmDouble>(lambdashy);
     3493   xDelete<IssmDouble>(n);
    34873494        xDelete<int>(doflist);
    34883495}/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.