Changeset 22511


Ignore:
Timestamp:
03/09/18 13:34:41 (7 years ago)
Author:
Mathieu Morlighem
Message:

CHG: removing some TripleMulitply

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

Legend:

Unmodified
Added
Removed

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

    r22485 r22511  
    569569        IssmDouble*    basis    = xNew<IssmDouble>(numnodes);
    570570        IssmDouble*    dbasis   = xNew<IssmDouble>(3*numnodes);
    571         IssmDouble*    B        = xNew<IssmDouble>(3*numnodes);
    572571        IssmDouble*    Bprime   = xNew<IssmDouble>(3*numnodes);
    573         IssmDouble     D[3][3]  = {0.};
    574572        IssmDouble     K[3][3];
    575573
    … …  
    600598
    601599                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     600                element->NodalFunctions(basis,gauss);
     601                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     602
    602603                D_scalar=gauss->weight*Jdet;
    603604                if(dt!=0.) D_scalar=D_scalar*dt;
    604605
    605606                /*Conduction: */
    606                 GetBConduct(B,element,xyz_list,gauss);
    607                 D[0][0]=D_scalar*kappa/rho_ice;
    608                 D[1][1]=D_scalar*kappa/rho_ice;
    609                 D[2][2]=D_scalar*kappa/rho_ice;
    610                 TripleMultiply(B,3,numnodes,1,
    611                                         &D[0][0],3,3,0,
    612                                         B,3,numnodes,0,
    613                                         &Ke->values[0],1);
     607                for(int i=0;i<numnodes;i++){
     608                        for(int j=0;j<numnodes;j++){
     609                                Ke->values[i*numnodes+j] += D_scalar*kappa/rho_ice*(
     610                                                        dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i]
     611                                                        );
     612                        }
     613                }
    614614
    615615                /*Advection: */
    616                 GetBAdvec(B,element,xyz_list,gauss);
    617                 GetBAdvecprime(Bprime,element,xyz_list,gauss);
    618616                vx_input->GetInputValue(&u,gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um;
    619617                vy_input->GetInputValue(&v,gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm;
    620618                vz_input->GetInputValue(&w,gauss); vzm_input->GetInputValue(&wm,gauss); vz=w-wm;
    621                 D[0][0]=D_scalar*vx;
    622                 D[1][1]=D_scalar*vy;
    623                 D[2][2]=D_scalar*vz;
    624                 TripleMultiply(B,3,numnodes,1,
    625                                         &D[0][0],3,3,0,
    626                                         Bprime,3,numnodes,0,
    627                                         &Ke->values[0],1);
     619                for(int i=0;i<numnodes;i++){
     620                        for(int j=0;j<numnodes;j++){
     621                                Ke->values[i*numnodes+j] += D_scalar*(
     622                                                        vx*dbasis[0*numnodes+j]*basis[i] + vy*dbasis[1*numnodes+j]*basis[i] +vz*dbasis[2*numnodes+j]*basis[i]
     623                                                        );
     624                        }
     625                }
    628626
    629627                /*Transient: */
    630628                if(dt!=0.){
    631629                        D_scalar=gauss->weight*Jdet;
    632                         element->NodalFunctions(basis,gauss);
    633                         TripleMultiply(basis,numnodes,1,0,
    634                                                 &D_scalar,1,1,0,
    635                                                 basis,1,numnodes,0,
    636                                                 &Ke->values[0],1);
     630                        for(int i=0;i<numnodes;i++){
     631                                for(int j=0;j<numnodes;j++){
     632                                        Ke->values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
     633                                }
     634                        }
    637635                        D_scalar=D_scalar*dt;
    638636                }
    … …  
    678676        xDelete<IssmDouble>(basis);
    679677        xDelete<IssmDouble>(dbasis);
    680         xDelete<IssmDouble>(B);
    681678        xDelete<IssmDouble>(Bprime);
    682679        delete gauss;

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

    r22506 r22511  
    242242
    243243        /*Intermediaries*/
    244         IssmDouble  D,Jdet;
     244        IssmDouble  Jdet;
    245245        IssmDouble *xyz_list = NULL;
    246246
    … …  
    250250        /*Initialize Element matrix and vectors*/
    251251        ElementMatrix* Ke     = element->NewElementMatrix(NoneApproximationEnum);
    252         IssmDouble*    B      = xNew<IssmDouble>(numnodes);
    253         IssmDouble*    Bprime = xNew<IssmDouble>(numnodes);
     252        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     253        IssmDouble*    dbasis = xNew<IssmDouble>(3*numnodes);
    254254
    255255        /*Retrieve all inputs and parameters*/
    … …  
    262262
    263263                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    264                 this->GetB(B,element,xyz_list,gauss);
    265                 this->GetBprime(Bprime,element,xyz_list,gauss);
    266                 D=gauss->weight*Jdet;
    267 
    268                 TripleMultiply(B,1,numnodes,1,
    269                                         &D,1,1,0,
    270                                         Bprime,1,numnodes,0,
    271                                         &Ke->values[0],1);
     264                element->NodalFunctions(basis,gauss);
     265                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     266
     267                for(int i=0;i<numnodes;i++){
     268                        for(int j=0;j<numnodes;j++){
     269                                Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
     270                                                        basis[j]*dbasis[2*numnodes+i]
     271                                                        );
     272                        }
     273                }
    272274        }
    273275
    … …  
    275277        delete gauss;
    276278        xDelete<IssmDouble>(xyz_list);
    277         xDelete<IssmDouble>(Bprime);
    278         xDelete<IssmDouble>(B);
     279        xDelete<IssmDouble>(dbasis);
     280        xDelete<IssmDouble>(basis);
    279281        return Ke;
    280282

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

    r22485 r22511  
    368368        int         stabilization;
    369369        IssmDouble  Jdet,dt,u,v,w,um,vm,wm,vel;
    370         IssmDouble  h,hx,hy,hz,vx,vy,vz;
     370        IssmDouble  h,hx,hy,hz,vx,vy,vz,D_scalar;
    371371        IssmDouble  tau_parameter,diameter;
    372         IssmDouble  D_scalar;
    373372        IssmDouble* xyz_list = NULL;
    374373
    … …  
    380379        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    381380        IssmDouble*    dbasis = xNew<IssmDouble>(3*numnodes);
    382         IssmDouble*    B      = xNew<IssmDouble>(3*numnodes);
    383381        IssmDouble*    Bprime = xNew<IssmDouble>(3*numnodes);
    384         IssmDouble     D[3][3]={0.};
    385382        IssmDouble     K[3][3];
    386383
    … …  
    409406
    410407                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     408                element->NodalFunctions(basis,gauss);
     409                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     410
    411411                D_scalar=gauss->weight*Jdet;
    412412                if(dt!=0.) D_scalar=D_scalar*dt;
    413413
    414414                /*Conduction: */
    415                 GetBConduct(B,element,xyz_list,gauss);
    416                 D[0][0]=D_scalar*kappa;
    417                 D[1][1]=D_scalar*kappa;
    418                 D[2][2]=D_scalar*kappa;
    419                 TripleMultiply(B,3,numnodes,1,
    420                                         &D[0][0],3,3,0,
    421                                         B,3,numnodes,0,
    422                                         &Ke->values[0],1);
     415                for(int i=0;i<numnodes;i++){
     416                        for(int j=0;j<numnodes;j++){
     417                                Ke->values[i*numnodes+j] += D_scalar*kappa*(
     418                                                        dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i]
     419                                                        );
     420                        }
     421                }
    423422
    424423                /*Advection: */
    425                 GetBAdvec(B,element,xyz_list,gauss);
    426                 GetBAdvecprime(Bprime,element,xyz_list,gauss);
    427424                vx_input->GetInputValue(&u,gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um;
    428425                vy_input->GetInputValue(&v,gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm;
    429426                vz_input->GetInputValue(&w,gauss); vzm_input->GetInputValue(&wm,gauss); vz=w-wm;
    430                 D[0][0]=D_scalar*vx;
    431                 D[1][1]=D_scalar*vy;
    432                 D[2][2]=D_scalar*vz;
    433                 TripleMultiply(B,3,numnodes,1,
    434                                         &D[0][0],3,3,0,
    435                                         Bprime,3,numnodes,0,
    436                                         &Ke->values[0],1);
     427                for(int i=0;i<numnodes;i++){
     428                        for(int j=0;j<numnodes;j++){
     429                                Ke->values[i*numnodes+j] += D_scalar*(
     430                                                        vx*dbasis[0*numnodes+j]*basis[i] + vy*dbasis[1*numnodes+j]*basis[i] +vz*dbasis[2*numnodes+j]*basis[i]
     431                                                        );
     432                        }
     433                }
    437434
    438435                /*Transient: */
    439436                if(dt!=0.){
    440437                        D_scalar=gauss->weight*Jdet;
    441                         element->NodalFunctions(basis,gauss);
    442                         TripleMultiply(basis,numnodes,1,0,
    443                                                 &D_scalar,1,1,0,
    444                                                 basis,1,numnodes,0,
    445                                                 &Ke->values[0],1);
     438                        for(int i=0;i<numnodes;i++){
     439                                for(int j=0;j<numnodes;j++){
     440                                        Ke->values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
     441                                }
     442                        }
    446443                        D_scalar=D_scalar*dt;
    447444                }
    … …  
    456453                        K[2][0]=h/(2.*vel)*fabs(vz*vx);  K[2][1]=h/(2.*vel)*fabs(vz*vy); K[2][2]=h/(2.*vel)*fabs(vz*vz);
    457454                        for(int i=0;i<3;i++) for(int j=0;j<3;j++) K[i][j] = D_scalar*K[i][j];
    458 
    459455                        GetBAdvecprime(Bprime,element,xyz_list,gauss);
    460456
    … …  
    486482        /*Clean up and return*/
    487483        xDelete<IssmDouble>(xyz_list);
    488         xDelete<IssmDouble>(B);
    489484        xDelete<IssmDouble>(Bprime);
    490485        xDelete<IssmDouble>(basis);
Note: See TracChangeset for help on using the changeset viewer.