Changeset 24429


Ignore:
Timestamp:
12/04/19 14:52:35 (5 years ago)
Author:
Mathieu Morlighem
Message:

CHG: trying to get rid of GetB

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

Legend:

Unmodified
Added
Removed

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

    r24335 r24429  
    116116
    117117        /*Initialize Element vector and other vectors*/
    118         ElementMatrix* Ke     = workelement->NewElementMatrix();
    119         IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
    120         IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
    121         IssmDouble*             dlsf   = xNew<IssmDouble>(dim);
    122         IssmDouble*             normal = xNew<IssmDouble>(dim);
    123    IssmDouble*          D           = xNewZeroInit<IssmDouble>(dim*dim);
     118        ElementMatrix *Ke     = workelement->NewElementMatrix();
     119        IssmDouble    *basis  = xNew<IssmDouble>(numnodes);
     120        IssmDouble    *dbasis = xNew<IssmDouble>(dim*numnodes);
     121        IssmDouble     dlsf[3];
     122        IssmDouble     normal[3];
    124123
    125124        /*Retrieve all inputs and parameters*/
    … …  
    143142
    144143                workelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
    145                 GetB(B,workelement,xyz_list,gauss,dim);
    146                 GetBprime(Bprime,workelement,xyz_list,gauss,dim);
     144                workelement->NodalFunctions(basis,gauss);
     145                workelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    147146
    148147                D_scalar=gauss->weight*Jdet;
    149148
    150149                if(extrapolatebydiffusion){
    151 
    152                         /* diffuse values outward only along the xy-plane*/
    153          for(int i=0;i<2;i++) D[i*dim+i] = D_scalar;
    154 
    155                         TripleMultiply(Bprime,dim,numnodes,1,
    156                                         D,dim,dim,0,
    157                                         Bprime,dim,numnodes,0,
    158                                         &Ke->values[0],1);
     150                        for(int i=0;i<numnodes;i++){
     151                                for(int j=0;j<numnodes;j++){
     152                                        Ke->values[i*numnodes+j] += D_scalar*(dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]);
     153                                }
     154                        }
    159155                }
    160156                else{
    … …  
    175171                                for(i=0;i<dim;i++)      normal[i]=0.;
    176172
    177                         for(row=0;row<dim;row++)
    178                                 for(col=0;col<dim;col++)
    179                                         if(row==col)
    180                                                 D[row*dim+col]=D_scalar*normal[row];
    181                                         else
    182                                                 D[row*dim+col]=0.;
    183                         TripleMultiply(B,dim,numnodes,1,
    184                                                 D,dim,dim,0,
    185                                                 Bprime,dim,numnodes,0,
    186                                                 &Ke->values[0],1);
     173                        for(int i=0;i<numnodes;i++){
     174                                for(int j=0;j<numnodes;j++){
     175                                        Ke->values[i*numnodes+j] += D_scalar*(normal[0]*dbasis[0*numnodes+j]*basis[i] + normal[1]*dbasis[1*numnodes+j]*basis[i]);
     176                                }
     177                        }
    187178
    188179                        /* stabilization */
    … …  
    195186                                h=sqrt(pow(hx*normal[0],2) + pow(hy*normal[1],2));
    196187                                kappa=h/2.+1.e-14;
    197                                 for(row=0;row<dim;row++)
    198                                         for(col=0;col<dim;col++)
    199                                                 if(row==col)
    200                                                         D[row*dim+col]=D_scalar*kappa;
    201                                                 else
    202                                                         D[row*dim+col]=0.;
    203                                 TripleMultiply(Bprime,dim,numnodes,1,
    204                                                         D,dim,dim,0,
    205                                                         Bprime,dim,numnodes,0,
    206                                                         &Ke->values[0],1);
     188
     189                                for(int i=0;i<numnodes;i++){
     190                                        for(int j=0;j<numnodes;j++){
     191                                                Ke->values[i*numnodes+j] += D_scalar*kappa*(dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]);
     192                                        }
     193                                }
    207194                        }
    208195                }
    … …  
    211198        /*Clean up and return*/
    212199        xDelete<IssmDouble>(xyz_list);
    213         xDelete<IssmDouble>(B);
    214         xDelete<IssmDouble>(Bprime);
    215         xDelete<IssmDouble>(D);
    216         xDelete<IssmDouble>(dlsf);
    217         xDelete<IssmDouble>(normal);
     200        xDelete<IssmDouble>(basis);
     201        xDelete<IssmDouble>(dbasis);
    218202        delete gauss;
    219203        if(extrapolationcase==0){workelement->DeleteMaterials(); delete workelement;};

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

    r24373 r24429  
    521521                                        basalelement->ElementSizes(&hx,&hy,&hz);
    522522                                        h=sqrt( pow(hx*w[0]/vel,2) + pow(hy*w[1]/vel,2) );
    523                                         IssmDouble D[9];
    524                                         for(row=0;row<dim;row++)
    525                                          for(col=0;col<dim;col++)
    526                                           D[row*dim+col] = D_scalar*h/(2.*vel)*w[row]*w[col];
    527                                         GetBprime(Bprime,basalelement,xyz_list,gauss);
    528 
    529                                         TripleMultiply(Bprime,dim,numnodes,1,
    530                                                                 &D[0],dim,dim,0,
    531                                                                 Bprime,dim,numnodes,0,
    532                                                                 &Ke->values[0],1);
     523                                        for(int i=0;i<numnodes;i++){
     524                                                for(int j=0;j<numnodes;j++){
     525                                                        Ke->values[i*numnodes+j] += D_scalar*h/(2.*vel)*w[i]*(w[j]*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + w[j]*dbasis[1*numnodes+j]*dbasis[1*numnodes+i]);
     526                                                }
     527                                        }
    533528                                  }
    534529                                break;
    … …  
    595590
    596591        return pe;
    597 }/*}}}*/
    598 void           LevelsetAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    599         /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
    600          * For node i, Bi can be expressed in the actual coordinate system
    601          * by:
    602          *       Bi=[ N ]
    603          *          [ N ]
    604          * where N is the finiteelement function for node i.
    605          *
    606          * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
    607          */
    608 
    609         /*Fetch number of nodes for this finite element*/
    610         int numnodes = element->GetNumberOfNodes();
    611 
    612         /*Get nodal functions*/
    613         IssmDouble* basis=xNew<IssmDouble>(numnodes);
    614         element->NodalFunctions(basis,gauss);
    615 
    616         /*Build B: */
    617         for(int i=0;i<numnodes;i++){
    618                 B[numnodes*0+i] = basis[i];
    619                 B[numnodes*1+i] = basis[i];
    620         }
    621 
    622         /*Clean-up*/
    623         xDelete<IssmDouble>(basis);
    624 }/*}}}*/
    625 void           LevelsetAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    626         /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
    627          * For node i, Bi' can be expressed in the actual coordinate system
    628          * by:
    629          *       Bi_prime=[ dN/dx ]
    630          *                [ dN/dy ]
    631          * where N is the finiteelement function for node i.
    632          *
    633          * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
    634          */
    635 
    636         /*Fetch number of nodes for this finite element*/
    637         int numnodes = element->GetNumberOfNodes();
    638 
    639         /*Get nodal functions derivatives*/
    640         IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
    641         element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    642 
    643         /*Build B': */
    644         for(int i=0;i<numnodes;i++){
    645                 Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
    646                 Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
    647         }
    648 
    649         /*Clean-up*/
    650         xDelete<IssmDouble>(dbasis);
    651 
    652592}/*}}}*/
    653593IssmDouble     LevelsetAnalysis::GetDistanceToStraight(IssmDouble* q, IssmDouble* s0, IssmDouble* s1){/*{{{*/

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

    r24335 r24429  
    2626                ElementMatrix* CreateKMatrix(Element* element);
    2727                ElementVector* CreatePVector(Element* element);
    28                 void           GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    29                 void           GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    3028                IssmDouble     GetDistanceToStraight(IssmDouble* q, IssmDouble* s0, IssmDouble* s1);
    3129                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);

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

    r24374 r24429  
    425425                if(stabilization==2){
    426426                        /*Streamline upwind*/
     427                        _assert_(dim==2);
    427428                        for(int i=0;i<numnodes;i++){
    428429                                for(int j=0;j<numnodes;j++){
    … …  
    499500        ElementMatrix* Ke     = element->NewElementMatrix();
    500501        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    501         IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
    502         IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
    503         IssmDouble     D[2][2];
     502        IssmDouble*    dbasis = xNew<IssmDouble>(3*numnodes);
    504503
    505504        /*Retrieve all inputs and parameters*/
    … …  
    517516                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    518517                element->NodalFunctions(basis,gauss);
     518                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    519519
    520520                vxaverage_input->GetInputValue(&vx,gauss);
    … …  
    522522
    523523                D_scalar=gauss->weight*Jdet;
    524                 TripleMultiply(basis,1,numnodes,1,
    525                                         &D_scalar,1,1,0,
    526                                         basis,1,numnodes,0,
    527                                         &Ke->values[0],1);
    528 
    529                 /*WARNING: B and Bprime are inverted compared to CG*/
    530                 GetB(Bprime,element,2,xyz_list,gauss);
    531                 GetBprime(B,element,2,xyz_list,gauss);
    532 
     524                for(int i=0;i<numnodes;i++){
     525                        for(int j=0;j<numnodes;j++){
     526                                Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j];
     527                        }
     528                }
     529
     530                /*WARNING: basis and dbasis are inverted compared to CG*/
    533531                D_scalar = - dt*gauss->weight*Jdet;
    534                 D[0][0]  = D_scalar*vx;
    535                 D[0][1]  = 0.;
    536                 D[1][0]  = 0.;
    537                 D[1][1]  = D_scalar*vy;
    538                 TripleMultiply(B,2,numnodes,1,
    539                                         &D[0][0],2,2,0,
    540                                         Bprime,2,numnodes,0,
    541                                         &Ke->values[0],1);
    542 
     532                for(int i=0;i<numnodes;i++){
     533                        for(int j=0;j<numnodes;j++){
     534                                Ke->values[i*numnodes+j] += D_scalar*(vx*dbasis[0*numnodes+i]*basis[j] + vy*dbasis[1*numnodes+i]*basis[j]);
     535                        }
     536                }
    543537        }
    544538
    … …  
    546540        xDelete<IssmDouble>(xyz_list);
    547541        xDelete<IssmDouble>(basis);
    548         xDelete<IssmDouble>(B);
    549         xDelete<IssmDouble>(Bprime);
     542        xDelete<IssmDouble>(dbasis);
    550543        delete gauss;
    551544        return Ke;
    … …  
    989982
    990983        /*Intermediaries */
    991         IssmDouble Jdet;
     984        IssmDouble Jdet,D_scalar;
    992985        IssmDouble vx,vy;
    993986        IssmDouble* xyz_list = NULL;
    … …  
    999992        /*Initialize Element vector and other vectors*/
    1000993        ElementMatrix* Ke     = element->NewElementMatrix();
    1001         IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
    1002         IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
     994        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     995        IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
    1003996        IssmDouble*    D      = xNewZeroInit<IssmDouble>(dim*dim);
    1004997
    … …  
    10141007
    10151008                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    1016                 GetB(B,element,dim,xyz_list,gauss);
    1017                 GetBprime(Bprime,element,dim,xyz_list,gauss);
     1009                element->NodalFunctions(basis,gauss);
     1010                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    10181011                vxaverage_input->GetInputValue(&vx,gauss);
    10191012                vyaverage_input->GetInputValue(&vy,gauss);
    10201013
    1021                 D[0*dim+0] = -gauss->weight*vx*Jdet;
    1022                 D[1*dim+1] = -gauss->weight*vy*Jdet;
    1023 
    1024                 TripleMultiply(B,dim,numnodes,1,
    1025                                         D,dim,dim,0,
    1026                                         Bprime,dim,numnodes,0,
    1027                                         &Ke->values[0],1);
    1028 
     1014                D_scalar = gauss->weight*Jdet;
     1015                for(int i=0;i<numnodes;i++){
     1016                        for(int j=0;j<numnodes;j++){
     1017                                Ke->values[i*numnodes+j] += -D_scalar*(vx*dbasis[0*numnodes+j]*basis[i] + vy*dbasis[1*numnodes+j]*basis[i]);
     1018                        }
     1019                }
    10291020        }
    10301021
    10311022        /*Clean up and return*/
    10321023        xDelete<IssmDouble>(xyz_list);
    1033         xDelete<IssmDouble>(B);
    1034         xDelete<IssmDouble>(Bprime);
    1035         xDelete<IssmDouble>(D);
     1024        xDelete<IssmDouble>(basis);
     1025        xDelete<IssmDouble>(dbasis);
    10361026        delete gauss;
    10371027        return Ke;

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

    r24335 r24429  
    124124
    125125        /*Initialize Element matrix and vectors*/
    126         ElementVector* pe    = element->NewElementVector();
    127         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    128         IssmDouble*    dvx   = xNew<IssmDouble>(dim);
    129         IssmDouble*    dvy   = xNew<IssmDouble>(dim);
    130         IssmDouble*    dvz   = xNew<IssmDouble>(dim);
     126        ElementVector *pe    = element->NewElementVector();
     127        IssmDouble    *basis = xNew<IssmDouble>(numnodes);
     128        IssmDouble     dvx[3];
     129        IssmDouble     dvy[3];
     130        IssmDouble     dvz[3];
    131131
    132132        Input2* vx_input=element->GetInput2(VxEnum);     _assert_(vx_input);
    … …  
    159159        xDelete<IssmDouble>(xyz_list);
    160160        xDelete<IssmDouble>(basis);
    161         xDelete<IssmDouble>(dvx);
    162         xDelete<IssmDouble>(dvy);
    163         xDelete<IssmDouble>(dvz);
    164161        return pe;
    165162}/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.