source: issm/oecreview/Archive/24684-25833/ISSM-25262-25263.diff

Last change on this file was 25834, checked in by Mathieu Morlighem, 4 years ago

CHG: added 24684-25833
File size: 10.7 KB

  • ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h

     
    2929                ElementMatrix* CreateJacobianMatrix(Element* element);
    3030                ElementMatrix* CreateKMatrix(Element* element);
    3131                ElementVector* CreatePVector(Element* element);
    32                 void           GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    33                 void           GetBprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    3432                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    3533                void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
    3634                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);

  • ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp

     
    453453        /*Initialize Element vector*/
    454454        ElementMatrix* Ke     = element->NewElementMatrix();
    455455        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    456         IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
    457         IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
     456        IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
    458457        IssmDouble*    D      = xNewZeroInit<IssmDouble>(dim*dim);
    459458
    460459        /*Retrieve all inputs and parameters*/
    461460        element->GetVerticesCoordinates(&xyz_list);
    462461        element->FindParam(&dt,TimesteppingTimeStepEnum);
    463         //printf("dt %f\n", dt);
    464462        element->FindParam(&stabilization,DamageStabilizationEnum);
    465463        Input2* vx_input = element->GetInput2(VxEnum); _assert_(vx_input);
    466464        Input2* vy_input = element->GetInput2(VyEnum); _assert_(vy_input);
    … …  
    478476
    479477                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    480478                element->NodalFunctions(basis,gauss);
     479                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    481480
    482481                vx_input->GetInputValue(&vx,gauss);
    483482                vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
    … …  
    489488                        vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
    490489                }
    491490
     491                /*Transient term*/
    492492                D_scalar=gauss->weight*Jdet;
    493                 TripleMultiply(basis,1,numnodes,1,
    494                                         &D_scalar,1,1,0,
    495                                         basis,1,numnodes,0,
    496                                         &Ke->values[0],1);
     493                for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j];
    497494
    498                 GetB(B,element,dim,xyz_list,gauss);
    499                 GetBprime(Bprime,element,dim,xyz_list,gauss);
    500 
    501495                dvxdx=dvx[0];
    502496                dvydy=dvy[1];
    503                 if(dim==3) dvzdz=dvz[2];
    504497                D_scalar=dt*gauss->weight*Jdet;
     498                if(dim==2){
     499                        for(int i=0;i<numnodes;i++){
     500                                for(int j=0;j<numnodes;j++){
     501                                        /*\phi_i \phi_j \nabla\cdot v*/
     502                                        Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j]*(dvxdx+dvydy);
     503                                        /*\phi_i v\cdot\nabla\phi_j*/
     504                                        Ke->values[i*numnodes+j] += D_scalar*basis[i]*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j]);
     505                                }
     506                        }
     507                }
     508                else{/*3D*/
     509                        _assert_(dim==3);
     510                        dvzdz=dvz[2];
     511                        for(int i=0;i<numnodes;i++){
     512                                for(int j=0;j<numnodes;j++){
     513                                        /*\phi_i \phi_j \nabla\cdot v*/
     514                                        Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j]*(dvxdx+dvydy+dvzdz);
     515                                        /*\phi_i v\cdot\nabla\phi_j*/
     516                                        Ke->values[i*numnodes+j] += D_scalar*basis[i]*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j] + vz*dbasis[2*numnodes+j]);
     517                                }
     518                        }
     519                }
    505520
    506                 D[0*dim+0]=D_scalar*dvxdx;
    507                 D[1*dim+1]=D_scalar*dvydy;
    508                 if(dim==3) D[2*dim+2]=D_scalar*dvzdz;
    509 
    510                 TripleMultiply(B,dim,numnodes,1,
    511                                         D,dim,dim,0,
    512                                         B,dim,numnodes,0,
    513                                         &Ke->values[0],1);
    514 
    515                 D[0*dim+0]=D_scalar*vx;
    516                 D[1*dim+1]=D_scalar*vy;
    517                 if(dim==3) D[2*dim+2]=D_scalar*vz;
    518 
    519                 TripleMultiply(B,dim,numnodes,1,
    520                                         D,dim,dim,0,
    521                                         Bprime,dim,numnodes,0,
    522                                         &Ke->values[0],1);
    523 
    524521                if(stabilization==2){
    525522                        if(dim==3){
    526523                                vel=sqrt(vx*vx+vy*vy+vz*vz)+1.e-8;
    … …  
    565562                                D[1*dim+0]=D_scalar*D[1*dim+0];
    566563                                D[0*dim+1]=D_scalar*D[0*dim+1];
    567564                                D[1*dim+1]=D_scalar*D[1*dim+1];
     565                                for(int i=0;i<numnodes;i++){
     566                                        for(int j=0;j<numnodes;j++){
     567                                                Ke->values[i*numnodes+j] += (
     568                                                                        dbasis[0*numnodes+i] *(D[0*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j]) +
     569                                                                        dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j])
     570                                                                        );
     571                                        }
     572                                }
    568573                        }
    569574                        else if(dim==3){
    570575                                D[0*dim+0]=D_scalar*D[0*dim+0];
    … …  
    576581                                D[0*dim+2]=D_scalar*D[0*dim+2];
    577582                                D[1*dim+2]=D_scalar*D[1*dim+2];
    578583                                D[2*dim+2]=D_scalar*D[2*dim+2];
     584                                for(int i=0;i<numnodes;i++){
     585                                        for(int j=0;j<numnodes;j++){
     586                                                Ke->values[i*numnodes+j] += (
     587                                                                        dbasis[0*numnodes+i] *(D[0*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j] + D[0*dim+2]*dbasis[2*numnodes+j]) +
     588                                                                        dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j] + D[1*dim+2]*dbasis[2*numnodes+j]) +
     589                                                                        dbasis[2*numnodes+i] *(D[2*dim+0]*dbasis[0*numnodes+j] + D[2*dim+1]*dbasis[1*numnodes+j] + D[2*dim+2]*dbasis[2*numnodes+j])
     590                                                                        );
     591                                        }
     592                                }
    579593                        }
    580                         TripleMultiply(Bprime,dim,numnodes,1,
    581                                                 D,dim,dim,0,
    582                                                 Bprime,dim,numnodes,0,
    583                                                 &Ke->values[0],1);
    584594                }
    585 
    586595        }
    587596
    588597        /*Clean up and return*/
    589598        xDelete<IssmDouble>(xyz_list);
    590599        xDelete<IssmDouble>(basis);
    591         xDelete<IssmDouble>(B);
    592         xDelete<IssmDouble>(Bprime);
     600        xDelete<IssmDouble>(dbasis);
    593601        xDelete<IssmDouble>(D);
    594602        delete gauss;
    595603        return Ke;
    … …  
    665673        delete gauss;
    666674        return pe;
    667675}/*}}}*/
    668 void           DamageEvolutionAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    669         /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*2.
    670          * For node i, Bi can be expressed in the actual coordinate system
    671          * by:
    672          *       Bi=[ N ]
    673          *          [ N ]
    674          * where N is the finiteelement function for node i.
    675          *
    676          * We assume B_prog has been allocated already, of size: 2x(1*numnodes)
    677          */
    678 
    679         /*Fetch number of nodes for this finite element*/
    680         int numnodes = element->GetNumberOfNodes();
    681 
    682         /*Get nodal functions*/
    683         IssmDouble* basis=xNew<IssmDouble>(numnodes);
    684         element->NodalFunctions(basis,gauss);
    685 
    686         /*Build B: */
    687         for(int i=0;i<numnodes;i++){
    688                 for(int j=0;j<dim;j++){
    689                         B[numnodes*j+i] = basis[i];
    690                 }
    691         }
    692 
    693         /*Clean-up*/
    694         xDelete<IssmDouble>(basis);
    695 }/*}}}*/
    696 void           DamageEvolutionAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    697         /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*2.
    698          * For node i, Bi' can be expressed in the actual coordinate system
    699          * by:
    700          *       Bi_prime=[ dN/dx ]
    701          *                [ dN/dy ]
    702          * where N is the finiteelement function for node i.
    703          *
    704          * We assume B' has been allocated already, of size: 3x(2*numnodes)
    705          */
    706 
    707         /*Fetch number of nodes for this finite element*/
    708         int numnodes = element->GetNumberOfNodes();
    709 
    710         /*Get nodal functions derivatives*/
    711         IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
    712         element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    713 
    714         /*Build B': */
    715         for(int i=0;i<numnodes;i++){
    716                 for(int j=0;j<dim;j++){
    717                         Bprime[numnodes*j+i] = dbasis[j*numnodes+i];
    718                 }
    719         }
    720 
    721         /*Clean-up*/
    722         xDelete<IssmDouble>(dbasis);
    723 
    724 }/*}}}*/
    725676void           DamageEvolutionAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    726677        element->GetSolutionFromInputsOneDof(solution,DamageDbarEnum);
    727678}/*}}}*/
    … …  
    787738
    788739        /*Initialize Element vector and other vectors*/
    789740        ElementMatrix* Ke     = element->NewElementMatrix();
    790         IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
    791         IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
    792         IssmDouble*    D      = xNewZeroInit<IssmDouble>(dim*dim);
     741        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     742        IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
    793743
    794744        /*Retrieve all inputs and parameters*/
    795745        element->GetVerticesCoordinates(&xyz_list);
    … …  
    802752                gauss->GaussPoint(ig);
    803753
    804754                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    805                 GetB(B,element,dim,xyz_list,gauss);
    806                 GetBprime(Bprime,element,dim,xyz_list,gauss);
     755                element->NodalFunctions(basis,gauss);
     756                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     757
    807758                vxaverage_input->GetInputValue(&vx,gauss);
    808759                vyaverage_input->GetInputValue(&vy,gauss);
    809760
    810                 D[0*dim+0] = -gauss->weight*vx*Jdet;
    811                 D[1*dim+1] = -gauss->weight*vy*Jdet;
    812 
    813                 TripleMultiply(B,dim,numnodes,1,
    814                                         D,dim,dim,0,
    815                                         Bprime,dim,numnodes,0,
    816                                         &Ke->values[0],1);
    817 
     761                for(int i=0;i<numnodes;i++){
     762                        for(int j=0;j<numnodes;j++){
     763                                /*\phi_i v\cdot\nabla\phi_j*/
     764                                Ke->values[i*numnodes+j] += -gauss->weight*Jdet*basis[i]*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j]);
     765                        }
     766                }
    818767        }
    819768
    820769        /*Clean up and return*/
    821770        xDelete<IssmDouble>(xyz_list);
    822         xDelete<IssmDouble>(B);
    823         xDelete<IssmDouble>(Bprime);
    824         xDelete<IssmDouble>(D);
     771        xDelete<IssmDouble>(basis);
     772        xDelete<IssmDouble>(dbasis);
    825773        delete gauss;
    826774        return Ke;
    827775}/*}}}*/
    … …  
    853801                element->NodalFunctions(basis,gauss);
    854802
    855803                D=gauss->weight*Jdet;
    856                 TripleMultiply(basis,1,numnodes,1,
    857                                         &D,1,1,0,
    858                                         basis,1,numnodes,0,
    859                                         &Me->values[0],1);
     804                for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Me->values[i*numnodes+j] += D*basis[i]*basis[j];
    860805        }
    861806
    862807        /*Clean up and return*/

  • ../trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp

     
    107107                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
    108108                basalelement->NodalFunctions(basis,gauss);
    109109                D=gauss->weight*Jdet;
    110 
    111                 TripleMultiply(basis,1,numnodes,1,
    112                                         &D,1,1,0,
    113                                         basis,1,numnodes,0,
    114                                         &Ke->values[0],1);
     110                for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D*basis[i]*basis[j];
    115111        }
    116112
    117113        /*Clean up and return*/

  • ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp

     
    336336
    337337                /*Transient term*/
    338338                D_scalar=gauss->weight*Jdet;
    339                 for(int i=0;i<numnodes;i++){
    340                         for(int j=0;j<numnodes;j++){
    341                                 Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j];
    342                         }
    343                 }
     339                for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j];
    344340
    345341                /*Advection terms*/
    346342                vxaverage_input->GetInputValue(&vx,gauss);
Note: See TracBrowser for help on using the repository browser.