Changeset 25264


Ignore:
Timestamp:
07/10/20 21:19:59 (5 years ago)
Author:
Mathieu Morlighem
Message:

CHG: removing TripleMultiply

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

Legend:

Unmodified
Added
Removed

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

    r25118 r25264  
    169169        ElementMatrix* Ke     = basalelement->NewElementMatrix(NoneApproximationEnum);
    170170        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    171         IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
    172         IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
     171        IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
    173172        IssmDouble*    D      = xNew<IssmDouble>(dim*dim);
    174173
    … …  
    189188                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
    190189                basalelement->NodalFunctions(basis,gauss);
     190                basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    191191
    192192                vx_input->GetInputValue(&vx,gauss);
    193193                if(dim==2) vy_input->GetInputValue(&vy,gauss);
    194194
     195                /*Transient term*/
    195196                D_scalar=gauss->weight*Jdet;
    196                 TripleMultiply(basis,1,numnodes,1,
    197                                         &D_scalar,1,1,0,
    198                                         basis,1,numnodes,0,
    199                                         &Ke->values[0],1);
    200 
    201                 GetB(B,basalelement,dim,xyz_list,gauss);
    202                 GetBprime(Bprime,basalelement,dim,xyz_list,gauss);
    203 
     197                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];
     198
     199                /*Advection terms*/
    204200                D_scalar=dt*gauss->weight*Jdet;
    205201                for(int i=0;i<dim*dim;i++) D[i]=0.;
    206                 D[0] = D_scalar*vx;
    207                 if(dim==2) D[1*dim+1] = D_scalar*vy;
    208 
    209                 TripleMultiply(B,dim,numnodes,1,
    210                                         D,dim,dim,0,
    211                                         Bprime,dim,numnodes,0,
    212                                         &Ke->values[0],1);
     202
     203                if(dim==1){
     204                        /*\phi_i v\cdot\nabla\phi_j*/
     205                        for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[i]*vx*dbasis[0*numnodes+j];
     206                }
     207                else{
     208                        _assert_(dim==2);
     209                        for(int i=0;i<numnodes;i++){
     210                                for(int j=0;j<numnodes;j++){
     211                                        /*\phi_i v\cdot\nabla\phi_j*/
     212                                        Ke->values[i*numnodes+j] += D_scalar*basis[i]*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j]);
     213                                }
     214                        }
     215                }
    213216
    214217                if(stabilization==2){
    … …  
    241244                if(stabilization==1 || stabilization==2){
    242245                        for(int i=0;i<dim*dim;i++) D[i]=D_scalar*D[i];
    243                         TripleMultiply(Bprime,dim,numnodes,1,
    244                                                 D,dim,dim,0,
    245                                                 Bprime,dim,numnodes,0,
    246                                                 &Ke->values[0],1);
     246                        if(dim==2){
     247                                for(int i=0;i<numnodes;i++){
     248                                        for(int j=0;j<numnodes;j++){
     249                                                Ke->values[i*numnodes+j] += (
     250                                                                        dbasis[0*numnodes+i] *(D[0*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j]) +
     251                                                                        dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j])
     252                                                                        );
     253                                        }
     254                                }
     255                        }
     256                        else{
     257                                for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += dbasis[0*numnodes+i]*D[0]*dbasis[0*numnodes+j];
     258                        }
    247259                }
    248260        }
    … …  
    251263        xDelete<IssmDouble>(xyz_list);
    252264        xDelete<IssmDouble>(basis);
    253         xDelete<IssmDouble>(B);
    254         xDelete<IssmDouble>(Bprime);
     265        xDelete<IssmDouble>(dbasis);
    255266        xDelete<IssmDouble>(D);
    256267        delete gauss;
    … …  
    332343        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    333344        return pe;
    334 
    335 }/*}}}*/
    336 void           FreeSurfaceBaseAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    337         /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*2.
    338          * For node i, Bi can be expressed in the actual coordinate system
    339          * by:
    340          *       Bi=[ N ]
    341          *          [ N ]
    342          * where N is the finiteelement function for node i.
    343          *
    344          * We assume B_prog has been allocated already, of size: 2x(1*numnodes)
    345          */
    346 
    347         /*Fetch number of nodes for this finite element*/
    348         int numnodes = element->GetNumberOfNodes();
    349 
    350         /*Get nodal functions*/
    351         IssmDouble* basis=xNew<IssmDouble>(numnodes);
    352         element->NodalFunctions(basis,gauss);
    353 
    354         /*Build B: */
    355         for(int i=0;i<numnodes;i++){
    356                 for(int j=0;j<dim;j++){
    357                         B[numnodes*j+i] = basis[i];
    358                 }
    359         }
    360 
    361         /*Clean-up*/
    362         xDelete<IssmDouble>(basis);
    363 }/*}}}*/
    364 void           FreeSurfaceBaseAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    365         /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*2.
    366          * For node i, Bi' can be expressed in the actual coordinate system
    367          * by:
    368          *       Bi_prime=[ dN/dx ]
    369          *                [ dN/dy ]
    370          * where N is the finiteelement function for node i.
    371          *
    372          * We assume B' has been allocated already, of size: 3x(2*numnodes)
    373          */
    374 
    375         /*Fetch number of nodes for this finite element*/
    376         int numnodes = element->GetNumberOfNodes();
    377 
    378         /*Get nodal functions derivatives*/
    379         IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
    380         element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    381 
    382         /*Build B': */
    383         for(int i=0;i<numnodes;i++){
    384                 for(int j=0;j<dim;j++){
    385                         Bprime[numnodes*j+i] = dbasis[j*numnodes+i];
    386                 }
    387         }
    388 
    389         /*Clean-up*/
    390         xDelete<IssmDouble>(dbasis);
    391345
    392346}/*}}}*/

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

    r24335 r25264  
    2626                ElementMatrix* CreateKMatrix(Element* element);
    2727                ElementVector* CreatePVector(Element* element);
    28                 void           GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    29                 void           GetBprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    3028                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    3129                void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);

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

    r25118 r25264  
    145145        ElementMatrix* Ke     = topelement->NewElementMatrix(NoneApproximationEnum);
    146146        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    147         IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
    148         IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
     147        IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
    149148        IssmDouble*    D      = xNew<IssmDouble>(dim*dim);
    150149
    … …  
    165164                topelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
    166165                topelement->NodalFunctions(basis,gauss);
     166                topelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    167167
    168168                vx_input->GetInputValue(&vx,gauss);
    169169                if(dim==2) vy_input->GetInputValue(&vy,gauss);
    170170
     171                /*Transient term*/
    171172                D_scalar=gauss->weight*Jdet;
    172                 TripleMultiply(basis,1,numnodes,1,
    173                                         &D_scalar,1,1,0,
    174                                         basis,1,numnodes,0,
    175                                         &Ke->values[0],1);
    176 
    177                 GetB(B,topelement,dim,xyz_list,gauss);
    178                 GetBprime(Bprime,topelement,dim,xyz_list,gauss);
    179 
     173                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];
     174
     175                /*Advection terms*/
    180176                D_scalar=dt*gauss->weight*Jdet;
    181177                for(int i=0;i<dim*dim;i++) D[i]=0.;
    … …  
    183179                if(dim==2) D[1*dim+1]=D_scalar*vy;
    184180
    185                 TripleMultiply(B,dim,numnodes,1,
    186                                         D,dim,dim,0,
    187                                         Bprime,dim,numnodes,0,
    188                                         &Ke->values[0],1);
     181                if(dim==1){
     182                        /*\phi_i v\cdot\nabla\phi_j*/
     183                        for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[i]*vx*dbasis[0*numnodes+j];
     184                }
     185                else{
     186                        _assert_(dim==2);
     187                        for(int i=0;i<numnodes;i++){
     188                                for(int j=0;j<numnodes;j++){
     189                                        /*\phi_i v\cdot\nabla\phi_j*/
     190                                        Ke->values[i*numnodes+j] += D_scalar*basis[i]*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j]);
     191                                }
     192                        }
     193                }
    189194
    190195                if(stabilization==2){
    … …  
    218223                if(stabilization==1 || stabilization==2){
    219224                        for(int i=0;i<dim*dim;i++) D[i]=D_scalar*D[i];
    220                         TripleMultiply(Bprime,dim,numnodes,1,
    221                                                 D,dim,dim,0,
    222                                                 Bprime,dim,numnodes,0,
    223                                                 &Ke->values[0],1);
     225                        if(dim==2){
     226                                for(int i=0;i<numnodes;i++){
     227                                        for(int j=0;j<numnodes;j++){
     228                                                Ke->values[i*numnodes+j] += (
     229                                                                        dbasis[0*numnodes+i] *(D[0*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j]) +
     230                                                                        dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j])
     231                                                                        );
     232                                        }
     233                                }
     234                        }
     235                        else{
     236                                for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += dbasis[0*numnodes+i]*D[0]*dbasis[0*numnodes+j];
     237                        }
    224238                }
    225239        }
    … …  
    228242        xDelete<IssmDouble>(xyz_list);
    229243        xDelete<IssmDouble>(basis);
    230         xDelete<IssmDouble>(B);
    231         xDelete<IssmDouble>(Bprime);
     244        xDelete<IssmDouble>(dbasis);
    232245        xDelete<IssmDouble>(D);
    233246        delete gauss;
    … …  
    306319
    307320}/*}}}*/
    308 void           FreeSurfaceTopAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    309         /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*2.
    310          * For node i, Bi can be expressed in the actual coordinate system
    311          * by:
    312          *       Bi=[ N ]
    313          *          [ N ]
    314          * where N is the finiteelement function for node i.
    315          *
    316          * We assume B_prog has been allocated already, of size: 2x(1*numnodes)
    317          */
    318 
    319         /*Fetch number of nodes for this finite element*/
    320         int numnodes = element->GetNumberOfNodes();
    321 
    322         /*Get nodal functions*/
    323         IssmDouble* basis=xNew<IssmDouble>(numnodes);
    324         element->NodalFunctions(basis,gauss);
    325 
    326         /*Build B: */
    327         for(int i=0;i<numnodes;i++){
    328                 for(int j=0;j<dim;j++){
    329                         B[numnodes*j+i] = basis[i];
    330                 }
    331         }
    332 
    333         /*Clean-up*/
    334         xDelete<IssmDouble>(basis);
    335 }/*}}}*/
    336 void           FreeSurfaceTopAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    337         /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*2.
    338          * For node i, Bi' can be expressed in the actual coordinate system
    339          * by:
    340          *       Bi_prime=[ dN/dx ]
    341          *                [ dN/dy ]
    342          * where N is the finiteelement function for node i.
    343          *
    344          * We assume B' has been allocated already, of size: 3x(2*numnodes)
    345          */
    346 
    347         /*Fetch number of nodes for this finite element*/
    348         int numnodes = element->GetNumberOfNodes();
    349 
    350         /*Get nodal functions derivatives*/
    351         IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
    352         element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    353 
    354         /*Build B': */
    355         for(int i=0;i<numnodes;i++){
    356                 for(int j=0;j<dim;j++){
    357                         Bprime[numnodes*j+i] = dbasis[j*numnodes+i];
    358                 }
    359         }
    360 
    361         /*Clean-up*/
    362         xDelete<IssmDouble>(dbasis);
    363 
    364 }/*}}}*/
    365321void           FreeSurfaceTopAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    366322           _error_("not implemented yet");

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

    r24335 r25264  
    2626                ElementMatrix* CreateKMatrix(Element* element);
    2727                ElementVector* CreatePVector(Element* element);
    28                 void           GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    29                 void           GetBprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    3028                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    3129                void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
Note: See TracChangeset for help on using the changeset viewer.