Changeset 25227


Ignore:
Timestamp:
07/07/20 21:37:04 (5 years ago)
Author:
Mathieu Morlighem
Message:

CHG: removing TripleMultiply from Hydro DC

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

Legend:

Unmodified
Added
Removed

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

    r25024 r25227  
    206206        /*Initialize Element vector*/
    207207        ElementMatrix* Ke     = basalelement->NewElementMatrix();
    208         IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
    209208        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    210         IssmDouble     D[2][2]={0.};
     209        IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
    211210
    212211        /*Retrieve all inputs and parameters*/
    … …  
    222221                gauss           ->GaussPoint(ig);
    223222                basalelement    ->JacobianDeterminant(&Jdet,xyz_list,gauss);
     223                basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     224                basalelement->NodalFunctions(basis,gauss);
    224225
    225226                epl_transmitivity = EplTransmitivity(basalelement,gauss,epl_thick_input,epl_head_input,base_input);
    … …  
    230231                //D_scalar=gauss->weight*Jdet;
    231232                if(dt!=0.) D_scalar=D_scalar*dt;
    232                 D[0][0]=D_scalar;
    233                 D[1][1]=D_scalar;
    234                 GetB(B,basalelement,xyz_list,gauss);
    235                 TripleMultiply(B,2,numnodes,1,
    236                                         &D[0][0],2,2,0,
    237                                         B,2,numnodes,0,
    238                                         &Ke->values[0],1);
     233                for(int i=0;i<numnodes;i++){
     234                        for(int j=0;j<numnodes;j++){
     235                                Ke->values[i*numnodes+j] += D_scalar*(dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]);
     236                        }
     237                }
    239238
    240239                /*Transient*/
    … …  
    243242                        D_scalar=epl_storing*gauss->weight*Jdet;
    244243                        //D_scalar=(epl_storing/epl_transmitivity)*gauss->weight*Jdet;
    245                         TripleMultiply(basis,numnodes,1,0,
    246                                                 &D_scalar,1,1,0,
    247                                                 basis,1,numnodes,0,
    248                                                 &Ke->values[0],1);
     244                        for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
    249245
    250246                        /*Transfer EPL part*/
    … …  
    252248                        D_scalar=dt*transfer*gauss->weight*Jdet;
    253249                        //D_scalar=dt*(transfer/epl_transmitivity)*gauss->weight*Jdet;
    254                         TripleMultiply(basis,numnodes,1,0,
    255                                                                                  &D_scalar,1,1,0,
    256                                                                                  basis,1,numnodes,0,
    257                                                                                  &Ke->values[0],1);
     250                        for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
    258251                }
    259252        }
    … …  
    262255        xDelete<IssmDouble>(xyz_list);
    263256        xDelete<IssmDouble>(basis);
    264         xDelete<IssmDouble>(B);
     257        xDelete<IssmDouble>(dbasis);
    265258        delete gauss;
    266259        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    … …  
    405398        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    406399        return pe;
    407 }/*}}}*/
    408 void HydrologyDCEfficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    409         /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*2.
    410          * For node i, Bi can be expressed in the actual coordinate system
    411          * by:
    412          *       Bi=[ dN/dx ]
    413          *          [ dN/dy ]
    414          * where N is the finiteelement function for node i.
    415          *
    416          * We assume B has been allocated already, of size: 3x(2*numnodes)
    417          */
    418 
    419         /*Fetch number of nodes for this finite element*/
    420         int numnodes = element->GetNumberOfNodes();
    421 
    422         /*Get nodal functions derivatives*/
    423         IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
    424         element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    425 
    426         /*Build B: */
    427         for(int i=0;i<numnodes;i++){
    428                 B[numnodes*0+i] = dbasis[0*numnodes+i];
    429                 B[numnodes*1+i] = dbasis[1*numnodes+i];
    430         }
    431 
    432         /*Clean-up*/
    433         xDelete<IssmDouble>(dbasis);
    434400}/*}}}*/
    435401void HydrologyDCEfficientAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/

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

    r24385 r25227  
    2929                ElementMatrix* CreateKMatrix(Element* element);
    3030                ElementVector* CreatePVector(Element* element);
    31                 void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    3231                void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    3332                void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);

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

    r25024 r25227  
    228228        /*Initialize Element vector*/
    229229        ElementMatrix* Ke     = basalelement->NewElementMatrix();
    230         IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
    231230        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    232         IssmDouble     D[2][2]= {0.};
     231        IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
    233232
    234233        /*Retrieve all inputs and parameters*/
    … …  
    252251                gauss          -> GaussPoint(ig);
    253252                basalelement   -> JacobianDeterminant(&Jdet,xyz_list,gauss);
     253                basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     254                basalelement->NodalFunctions(basis,gauss);
    254255
    255256                sediment_transmitivity = SedimentTransmitivity(basalelement,gauss,sed_head_input,base_input,SedTrans_input);
    … …  
    260261                //D_scalar=gauss->weight*Jdet;
    261262                if(dt!=0.) D_scalar=D_scalar*dt;
    262                 D[0][0]=D_scalar;
    263                 D[1][1]=D_scalar;
    264                 GetB(B,basalelement,xyz_list,gauss);
    265                 TripleMultiply(B,2,numnodes,1,
    266                                                                          &D[0][0],2,2,0,
    267                                                                          B,2,numnodes,0,
    268                                                                          &Ke->values[0],1);
     263                for(int i=0;i<numnodes;i++){
     264                        for(int j=0;j<numnodes;j++){
     265                                Ke->values[i*numnodes+j] += D_scalar*(dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]);
     266                        }
     267                }
    269268
    270269                /*Transient*/
    271270                if(dt!=0.){
    272                         basalelement->NodalFunctions(&basis[0],gauss);
    273271                        D_scalar=sediment_storing*gauss->weight*Jdet;
    274272                        //D_scalar=(sediment_storing/sediment_transmitivity)*gauss->weight*Jdet;
    275                         TripleMultiply(basis,numnodes,1,0,
    276                                                                                  &D_scalar,1,1,0,
    277                                                                                  basis,1,numnodes,0,
    278                                                                                  &Ke->values[0],1);
     273                        for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
    279274
    280275                        /*Transfer EPL part*/
    … …  
    282277                                if(active_element){
    283278                                        transfer=GetHydrologyKMatrixTransfer(basalelement);
    284                                         basalelement->NodalFunctions(&basis[0],gauss);
    285279                                        D_scalar=dt*transfer*gauss->weight*Jdet;
    286280                                        //D_scalar=dt*(transfer/sediment_transmitivity)*gauss->weight*Jdet;
    287                                         TripleMultiply(basis,numnodes,1,0,
    288                                                                                                  &D_scalar,1,1,0,
    289                                                                                                  basis,1,numnodes,0,
    290                                                                                                  &Ke->values[0],1);
     281                                        for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
    291282                                }
    292283                        }
    … …  
    295286        /*Clean up and return*/
    296287        xDelete<IssmDouble>(xyz_list);
    297         xDelete<IssmDouble>(B);
    298288        xDelete<IssmDouble>(basis);
     289        xDelete<IssmDouble>(dbasis);
    299290        delete gauss;
    300291        if(domaintype!=Domain2DhorizontalEnum){
    … …  
    473464        }
    474465        return pe;
    475 }/*}}}*/
    476 void HydrologyDCInefficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    477         /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*2.
    478          * For node i, Bi can be expressed in the actual coordinate system
    479          * by:
    480          *       Bi=[ dN/dx ]
    481          *          [ dN/dy ]
    482          * where N is the finiteelement function for node i.
    483          *
    484          * We assume B has been allocated already, of size: 3x(2*numnodes)
    485          */
    486 
    487         /*Fetch number of nodes for this finite element*/
    488         int numnodes = element->GetNumberOfNodes();
    489 
    490         /*Get nodal functions derivatives*/
    491         IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
    492         element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    493 
    494         /*Build B: */
    495         for(int i=0;i<numnodes;i++){
    496                 B[numnodes*0+i] = dbasis[0*numnodes+i];
    497                 B[numnodes*1+i] = dbasis[1*numnodes+i];
    498         }
    499 
    500         /*Clean-up*/
    501         xDelete<IssmDouble>(dbasis);
    502466}/*}}}*/
    503467void HydrologyDCInefficientAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/

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

    r24335 r25227  
    3333
    3434                /*Intermediaries*/
    35                 void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    3635                IssmDouble SedimentStoring(Element* element, Gauss* gauss, Input2* sed_head_input, Input2* base_input);
    3736                IssmDouble SedimentTransmitivity(Element* element,Gauss* gauss,Input2* sed_head_input, Input2* base_input,Input2* SedTrans_input);
Note: See TracChangeset for help on using the changeset viewer.