/Users/jack/Code/basso_dev/inc/basic_fem_operations.h

Go to the documentation of this file.
00001 
00015 #ifndef _BASIC_FEM_OPERATIONS_BASSO_H_
00016 #define _BASIC_FEM_OPERATIONS_BASSO_H_
00017 
00018 // std includes
00019 #include <map>
00020 
00021 // Gmm++ includes
00022 #include "gmm.h"
00023 #include "gmm_superlu_interface.h"
00024 
00025 // Basso includes
00026 #include "basso.h"
00027 #include "basic_functions.h"
00028 #include "Node.h"
00029 
00030 using namespace std;
00031 using namespace gmm;
00032 
00033 namespace basso {
00034 
00040         template < class NODEARRAY >
00041         void form_nodemap( NODEARRAY &nodes, map<int,Node*> &nodemap )
00042         {
00043                 typename NODEARRAY::iterator nItr;      
00044                 for ( nItr=nodes.begin(); nItr!=nodes.end(); ++nItr )
00045                         nodemap[ nItr->Id() ] = &(*nItr);
00046         }
00047 
00053         void form_vect_shapefunct( const nArray &N, int sdim,  nMatrix &Nv )
00054         {
00055                 int nn=N.size();
00056 #ifdef ALLOW_DYNAMIC_RESIZE
00057                 if ( mat_nrows(Nv)!=nn || mat_ncols(Nv)!=sdim )
00058                         resize( Nv, nn, sdim );
00059 #elif EXPLICIT_BOUNDS_CHECK 
00060                 if ( mat_nrows(Nv)!=nn || mat_ncols(Nv)!=sdim )
00061                         ErrorMessage("form_vect_shapefunct","incompatible dimensions");
00062 #endif
00063                 clear(Nv);
00064                 for ( int i=0; i<nn; ++i )
00065                         for ( int s=0; s<sdim; ++s )
00066                                 Nv(i+s,s)=N[i];
00067         }
00068 
00073         void form_bmatrix( const nMatrix &DNa_x, nMatrix &B )
00074         {
00075                 int nn=mat_nrows(DNa_x), sdim=mat_ncols(DNa_x);
00076                 int vdim=1;
00077                 if ( sdim==2 )
00078                         vdim=3;
00079                 else if ( sdim==3 )
00080                         vdim=6;
00081 #ifdef ALLOW_DYNAMIC_RESIZE
00082                 if ( mat_nrows(B)!=vdim || mat_ncols(B)!=sdim*nn )
00083                         resize( B, vdim, sdim*nn );
00084 #elif EXPLICIT_BOUNDS_CHECK 
00085                 if ( mat_nrows(B)!=vdim || mat_ncols(B)!=sdim*nn )
00086                         ErrorMessage("form_vect_shapefunct","incompatible dimensions");
00087 #endif  
00088                 clear(B);
00089                 if ( sdim==1 )
00090                         copy(transposed(DNa_x),B);
00091                 else if ( sdim==2 ) 
00092                         for ( int i=0; i<nn; ++i ) {
00093                                 B(0,2*i)   = DNa_x(i,0);
00094                                 B(1,2*i+1) = DNa_x(i,1);
00095                                 B(2,2*i)   = DNa_x(i,1);
00096                                 B(2,2*i+1) = DNa_x(i,0);        
00097                         }
00098         
00099                 else // sdim==3 
00100                         for ( int i=0; i<nn; ++i ) {
00101                                 B(0,3*i)   = DNa_x(i,0);
00102                                 B(1,3*i+1) = DNa_x(i,1);
00103                                 B(2,3*i+2) = DNa_x(i,2);
00104                                 B(3,3*i+1) = DNa_x(i,2);
00105                                 B(3,3*i+2) = DNa_x(i,1);
00106                                 B(4,3*i+2) = DNa_x(i,0);
00107                                 B(4,3*i)   = DNa_x(i,2);
00108                                 B(5,3*i)   = DNa_x(i,1);
00109                                 B(5,3*i+1) = DNa_x(i,0);        
00110                         }
00111         }
00112 
00122         template < class NODEARRAY, class CONNVECT, class COORDMAT >
00123         void get_element_coordinates( const NODEARRAY &node, int sdim,
00124                                                                                 const CONNVECT &econn, COORDMAT &cmat )
00125         {
00126                 int nn=vect_size(econn);
00127 #ifdef ALLOW_DYNAMIC_RESIZE
00128                 if ( mat_nrows(cmat)!=nn || mat_ncols(cmat)!=sdim )
00129                         resize( cmat, nn, sdim );
00130 #elif EXPLICIT_BOUNDS_CHECK 
00131                 if ( mat_nrows(cmat)!=nn || mat_ncols(cmat)!=sdim )
00132                         ErrorMessage("get_element_coordinates","incompatible dimensions");
00133 #endif          
00134                 for ( int n=0; n<nn; ++n ) {
00135                         int nn = econn[n];
00136                         for ( int i=0; i<sdim; ++i )
00137                                 cmat(n,i) = node[nn].x(i);
00138                 }
00139         } 
00140 
00148         template< class OBJ1, class OBJ2, class INDX >
00149         void scatter_vector( const OBJ1 &fe, OBJ2 &F, const INDX &indx )
00150         {
00151 #ifdef EXPLICIT_BOUNDS_CHECK
00152                 if ( vect_size(fe)!=vect_size(indx)  )
00153                         ErrorMessage("scatter_vector","incompatible dimensions");
00154 #endif          
00155                 for ( int i=0; i<vect_size(indx); ++i )
00156                         F[ indx[i] ] += fe[i];
00157         }
00158         
00167         template< class OBJ1, class OBJ2, class INDX >
00168         void scatter_matrix( const OBJ1 &ke, OBJ2 &K, const INDX &ri, const INDX &ci )
00169         {
00170 #ifdef EXPLICIT_BOUNDS_CHECK
00171                 if ( mat_nrows(ke)!=vect_size(ri) || mat_ncols(ke)!=vect_size(ci) )
00172                         ErrorMessage("scatter_matrix","incompatible dimensions");
00173 #endif
00174                 for ( int i=0; i<vect_size(ri); ++i )
00175                         for ( int j=0; j<vect_size(ci); ++j ) 
00176                                 K(ri[i],ci[j])+=ke(i,j);
00177         }
00178         
00188         template< class INDX, class NVECT >
00189         void superlu_solve_constrained( DynamicCSRMatrix &K, nArray &x, const INDX &ifix, NVECT &ival )
00190         {
00191                 int nc=ifix.size(), ndof=mat_nrows(K);
00192 
00193                 // get reaction equations
00194                 DynamicCSRMatrix Kreac(nc,ndof);
00195                 copy( sub_matrix( K, unsorted_sub_index(ifix), sub_interval(0,ndof) ), Kreac );
00196 
00197                 // modify  rhs f := -K(:,ifix)*ival + f 
00198                 mult_add( sub_matrix( K, sub_interval(0,ndof), sub_index(ifix) ), scaled(ival,-1.0), x ); 
00199 
00200                 // zero out rows and columns
00201                 for ( int i=0; i<nc; ++i )
00202                 {
00203                         int ii=ifix[i];
00204                         for ( int jj=0; jj<ndof; ++jj )
00205                         {
00206                                 if ( ii != jj)
00207                                 {
00208                                         K(ii,jj)=0.0; 
00209                                         K(jj,ii)=0.0;
00210                                 }
00211                         }
00212                         x[ii] = ival[i]*K(ii,ii);
00213                 }
00214 
00215                 // solve the system
00216                 CSRMatrix Kcsr(ndof,ndof);
00217                 clean(K,1e-12);
00218                 copy(K,Kcsr);
00219                 Numeric condest;
00220                 SuperLU_solve( Kcsr, x, x, condest, 1 );
00221 
00222                 // get reaction forces
00223                 mult( Kreac, x, ival );
00224         }
00225 
00226         
00235         template< class INDX >
00236         void superlu_solve_constrained( DynamicCSRMatrix &K, nArray &x, const INDX &ifix )
00237         {
00238                 int nc=ifix.size(), ndof=mat_nrows(K);
00239 
00240                 // zero out rows and columns
00241                 typename INDX::const_iterator ifixItr;
00242                 ifixItr=ifix.begin();
00243                 for ( int i=0; i<nc; ++i, ++ifixItr )
00244                 {
00245                         int ii = *ifixItr;
00246                         for ( int jj=0; jj<ndof; ++jj )
00247                         {
00248                                 if ( ii != jj)
00249                                 {
00250                                         K(ii,jj)=0.0; 
00251                                         K(jj,ii)=0.0;
00252                                 }
00253                         }
00254                         x[ii] = 0.0;
00255                 }
00256 
00257                 // solve the system
00258                 CSRMatrix Kcsr(ndof,ndof);
00259                 clean(K,1e-12);
00260                 copy(K,Kcsr);
00261                 Numeric condest;
00262                 SuperLU_solve( Kcsr, x, x, condest, 1 );
00263 
00264         }
00265         
00266         /****************************************************************************************************** 
00267          ***                                                                                                ***
00268          ***                Functions to compute determinants and inverse of small matrices                 ***
00269          ***                                                                                                ***
00270          ******************************************************************************************************/
00271 
00274         Numeric det2by2( const nMatrix &A )
00275         {
00276                 return A(0,0)*A(1,1)-A(0,1)*A(1,0);
00277         }
00278         
00281         Numeric det3by3( const nMatrix &A )
00282         {
00283                 return A(0,0)*(A(1,1)*A(2,2)-A(1,2)*A(2,1))+A(0,1)*(A(1,2)*A(2,0)-A(1,0)*A(2,2))
00284                         +A(0,2)*(A(1,0)*A(2,1)-A(1,1)*A(2,0));
00285         }       
00286         
00290         Numeric inv2by2( nMatrix &A )
00291         {
00292                 Numeric d=det2by2(A), invd=1.0/d, A00=A(0,0);
00293                 A(0,0)=invd*A(1,1); A(0,1)=-invd*A(0,1);
00294                 A(1,0)=-invd*A(1,0); A(1,1)=invd*A00;
00295                 return d;
00296         }
00297         
00301         Numeric inv3by3( nMatrix &A )
00302         {
00303                 Numeric d=det3by3(A), invd=1.0/d;
00304                 nMatrix Ainv(3,3);
00305                 Ainv(0,0)=invd*( -A(1,2)*A(2,1) + A(1,1)*A(2,2) ); 
00306                 Ainv(0,1)=invd*( A(0,2)*A(2,1) - A(0,1)*A(2,2) ); 
00307                 Ainv(0,2)=invd*(-A(0,2)*A(1,1) + A(0,1)*A(1,2) );       
00308                 Ainv(1,0)=invd*( A(1,2)*A(2,0) - A(1,0)*A(2,2) );
00309                 Ainv(1,1)=invd*( -A(0,2)*A(2,0) + A(0,0)*A(2,2) ); 
00310                 Ainv(1,2)=invd*( A(0,2)*A(1,0) - A(0,0)*A(1,2) );               
00311                 Ainv(2,0)=invd*( -A(1,1)*A(2,0) + A(1,0)*A(2,1) ); 
00312                 Ainv(2,1)=invd*( A(0,1)*A(2,0) - A(0,0)*A(2,1) ); 
00313                 Ainv(2,2)=invd*( -A(0,1)*A(1,0) + A(0,0)*A(1,1) );
00314                 copy(Ainv,A);
00315                 return d;
00316         }
00317         
00318         
00319         /****************************************************************************************************** 
00320         ***                                                                                                ***
00321         ***                Functions to compute spacial gradients, B matrix and Jacobians                  ***
00322         ***                                                                                                ***
00323         ******************************************************************************************************/ 
00324         
00336         void jacobian_mat( const nMatrix &coord, const nMatrix &dNxi, nMatrix &jac )
00337         {
00338                 int nn=mat_nrows(dNxi), sdim=mat_ncols(dNxi);           
00339                 mult( transposed(sub_matrix(coord, sub_interval(0, nn), sub_interval(0, sdim))), dNxi, jac );
00340         }
00347         Numeric jacobian( const nMatrix &coord, const nMatrix &dNxi )
00348         {
00349                 int sdim=mat_ncols(dNxi);
00350                 nMatrix jac(sdim,sdim);
00351                 jacobian_mat(coord,dNxi,jac);
00352                 if ( sdim==1 )
00353                         return jac(0,0);
00354                 else if ( sdim==2 )
00355                         return det2by2(jac); 
00356                 else if ( sdim==3 )
00357                         return det3by3(jac);
00358                 else
00359                         return lu_det(jac);
00360         }       
00361         
00369         Numeric jacobian12( const nMatrix &coord, const nMatrix &dNxi )
00370         {
00371                 int nn=mat_nrows(dNxi), sdim=2, edim=1;
00372                         
00373                 nMatrix jac(sdim,sdim);
00374                 clear(jac);
00375                 for ( int i=0; i<nn; ++i ) {
00376                         jac(0,0) += coord(i,0)*dNxi(i,0);  // x,xi  
00377                         jac(1,0) += coord(i,1)*dNxi(i,0);  // y,xi
00378                 }
00379                 Numeric j=det2by2(jac);
00380                 return sqrt(j);  // sqrt( dxdxi ^2 + dydxi ^2 )
00381         }
00382         
00390         Numeric jacobian13( const nMatrix &coord, const nMatrix &dNxi )
00391         {
00392                 int nn=mat_nrows(dNxi), sdim=3, edim=1;
00393                 
00394                 nArray Exi(3);
00395                 clear(Exi);
00396                 for ( int i=0; i<nn; ++i ) {
00397                         Exi[0] += coord(i,0)*dNxi(i,0);  // x,xi  
00398                         Exi[1] += coord(i,1)*dNxi(i,0);  // y,xi  
00399                         Exi[2] += coord(i,2)*dNxi(i,0);  // z,xi
00400                 }               
00401 
00402                 return sqrt( pow(Exi[0],2) + pow(Exi[1],2) + pow(Exi[2],2) );
00403         }
00411         Numeric jacobian23( const nMatrix &coord, const nMatrix &dNxi )
00412         {
00413                 int nn=mat_nrows(dNxi), sdim=3, edim=2;
00414                 nArray Exi(3), Eeta(3), Ezeta(3);
00415                 clear(Exi); clear(Eeta);
00416                 for ( int i=0; i<nn; ++i ) {
00417                         Exi[0] += coord(i,0)*dNxi(i,0);  // x,xi  
00418                         Exi[1] += coord(i,1)*dNxi(i,0);  // y,xi  
00419                         Exi[2] += coord(i,2)*dNxi(i,0);  // z,xi
00420                         Eeta[0] += coord(i,0)*dNxi(i,1);  // x,eta  
00421                         Eeta[1] += coord(i,1)*dNxi(i,1);  // y,eta 
00422                         Eeta[2] += coord(i,2)*dNxi(i,1);  // z,eta
00423                 }
00424 
00425                 cross_prod(Exi,Eeta,Ezeta);
00426                 return vect_norm2( Ezeta );
00427         }
00439         Numeric grad_shape( const nMatrix &coord, const nMatrix &dNxi, nMatrix &dNx )
00440         {
00441                 int nn=mat_nrows(dNxi), sdim=mat_ncols(dNxi);
00442 #ifdef ALLOW_DYNAMIC_RESIZE
00443                 if ( mat_nrows(dNx)!=nn || mat_ncols(dNx)!=sdim  ) resize( dNx, nn, sdim );
00444 #elif EXPLICIT_BOUNDS_CHECK 
00445                 if ( mat_nrows(dNx)!=nn || mat_ncols(dNx)!=sdim ) ErrorMessage("grad_shape","dNx of incorrect dimension");
00446                 if ( mat_nrows(coord)!=nn || mat_ncols(coord)>=sdim ) ErrorMessage("grad_shape","coord of incorrect dimension");
00447 #endif          
00448                 nMatrix invj(sdim,sdim);
00449                 Numeric j;
00450                 
00451                 mult(transposed(coord),dNxi,invj);
00452                 
00453                 if ( sdim==1 ) {
00454                         j=invj(0,0); 
00455                         invj(0,0)=1.0/j;
00456                 }
00457                 else if ( sdim==2 )
00458                         j=inv2by2(invj); 
00459                 else if ( sdim==3 )
00460                         j=inv3by3(invj); 
00461                 else
00462                         j=lu_inverse(invj);
00463                         
00464                 mult(dNxi,invj,dNx);
00465                 return j;
00466         }
00467 
00481         Numeric grad_shape( const Array<Node> &nodes, const iArray &conn, const nMatrix &dNxi, nMatrix &invj, nMatrix &dNx )
00482         {
00483                 int nn=mat_nrows(dNxi), sdim=mat_ncols(dNxi);
00484                 
00485                 Numeric j;
00486                 
00487                 clear(invj);
00488                 for  ( int I=0; I<nn; ++I )
00489                         for ( int i=0; i<sdim; ++i )
00490                                 for ( int j=0; j<sdim; ++j )
00491                                         invj(i,j) += nodes[ conn[I] ].x(i) * dNxi(I,j);
00492 
00493                 if ( sdim==1 ) {
00494                         j=invj(0,0); 
00495                         invj(0,0)=1.0/j;
00496                 }
00497                 else if ( sdim==2 )
00498                         j=inv2by2(invj); 
00499                 else if ( sdim==3 )
00500                         j=inv3by3(invj); 
00501                 else
00502                         j=lu_inverse(invj);
00503 
00504                 mult(dNxi,invj,dNx);
00505                 return j;
00506         }
00507         
00514         Numeric grad_shape12( const nMatrix &coord, const nMatrix &dNxi, nMatrix &dNx )
00515         {
00516                 int nn=mat_nrows(dNxi), sdim=2, edim=1;
00517 #ifdef ALLOW_DYNAMIC_RESIZE
00518                 if ( mat_nrows(dNx)!=nn || mat_ncols(dNx)!=sdim  ) resize( dNx, nn, sdim );
00519 #elif EXPLICIT_BOUNDS_CHECK 
00520                 if ( mat_nrows(dNx)!=nn || mat_ncols(dNx)!=sdim ) ErrorMessage("gradshape12","dNx of incorrect dimension");
00521                 if ( mat_nrows(coord)!=nn || mat_ncols(coord)>=sdim ) ErrorMessage("gradshape12","coord of incorrect dimension");
00522 #endif                          
00523                 nMatrix invj(sdim,sdim);
00524                 clear(invj);
00525                 for ( int i=0; i<nn; ++i ) {
00526                         invj(0,0) += coord(i,0)*dNxi(i,0);  // x,xi  
00527                         invj(1,0) += coord(i,1)*dNxi(i,0);  // y,xi
00528                 }
00529                 invj(0,1) = -invj(0,0); // x,eta=-y,xi
00530                 invj(1,1) =  invj(1,0); // y,eta=x,xi
00531                 
00532                 Numeric j=inv2by2(invj);
00533                 for ( int i=0; i<nn; ++i)
00534                 {
00535                         dNx(i,0)=dNxi(i,0)*invj(0,0);  // dNIdxi dxidx
00536                         dNx(i,1)=dNxi(i,0)*invj(0,1);  // dNIdxi dxidy
00537                 }
00538                 
00539                 return sqrt(j);  // sqrt( dxdxi ^2 + dydxi ^2 )
00540         }
00541         
00542 
00549         Numeric grad_shape13( const nMatrix &coord, const nMatrix &dNxi, nMatrix &dNx )
00550         {
00551                 int nn=mat_nrows(dNxi), sdim=3, edim=1;
00552 #ifdef ALLOW_DYNAMIC_RESIZE
00553                 if ( mat_nrows(dNx)!=nn || mat_ncols(dNx)!=sdim  ) resize( dNx, nn, sdim );
00554 #elif EXPLICIT_BOUNDS_CHECK 
00555                 if ( mat_nrows(dNx)!=nn || mat_ncols(dNx)!=sdim ) ErrorMessage("gradshape13","dNx of incorrect dimension");
00556                 if ( mat_nrows(coord)!=nn || mat_ncols(coord)>=sdim ) ErrorMessage("gradshape13","coord of incorrect dimension");
00557 #endif          
00558                 nMatrix invj(sdim,sdim);
00559                 
00560                 nArray Exi(3), Eeta(3), Ezeta(3);
00561                 clear(Exi);
00562                 for ( int i=0; i<nn; ++i ) {
00563                         Exi[0] += coord(i,0)*dNxi(i,0);  // x,xi  
00564                         Exi[1] += coord(i,1)*dNxi(i,0);  // y,xi  
00565                         Exi[2] += coord(i,2)*dNxi(i,0);  // z,xi
00566                 }               
00567                 
00568                 if ( Exi[0]!=0.0 || Exi[1]!=0.0 ) {
00569                         Eeta[0]=Exi[1]; Eeta[1]=-Exi[0]; Eeta[2]=0.0;  // Eeta = Exi x e3
00570                 }
00571                 else  {// Exi parallel to e3
00572                         Eeta[0]=Exi[2]; Eeta[1]=0.0; Eeta[2]=-Exi[0];  // Eeta = e2 x Exi
00573                 }
00574                 cross_prod(Exi,Eeta,Ezeta);
00575                                 
00576                 invj(0,0)=Exi[0];  invj(0,1)=Eeta[0]; invj(0,2)=Ezeta[0];
00577                 invj(1,0)=Exi[1];  invj(1,1)=Eeta[1]; invj(1,2)=Ezeta[1];
00578                 invj(2,0)=Exi[2];  invj(2,1)=Eeta[2]; invj(2,2)=Ezeta[2];
00579                 
00580                 Numeric j=inv3by3(invj);
00581                 for ( int i=0; i<nn; ++i)
00582                 {
00583                         dNx(i,0)=dNxi(i,0)*invj(0,0); // dNIdxi dxidx
00584                         dNx(i,1)=dNxi(i,0)*invj(0,1); // dNIdxi dxidy
00585                         dNx(i,2)=dNxi(i,0)*invj(0,2); // dNIdxi dxidz
00586                 }
00587                 
00588                 return sqrt( pow(Exi[0],2) + pow(Exi[1],2) + pow(Exi[2],2) );
00589         }       
00590         
00591 
00598         Numeric grad_shape23( const nMatrix &coord, const nMatrix &dNxi, nMatrix &dNx )
00599         {
00600                 int nn=mat_nrows(dNxi), sdim=3, edim=2;
00601 #ifdef ALLOW_DYNAMIC_RESIZE
00602                 if ( mat_nrows(dNx)!=nn || mat_ncols(dNx)!=sdim  ) resize( dNx, nn, sdim );
00603 #elif EXPLICIT_BOUNDS_CHECK 
00604                 if ( mat_nrows(dNx)!=nn || mat_ncols(dNx)!=sdim ) ErrorMessage("gradshape23","dNx of incorrect dimension");
00605                 if ( mat_nrows(coord)!=nn || mat_ncols(coord)>=sdim ) ErrorMessage("gradshape23","coord of incorrect dimension");
00606 #endif          
00607                 nMatrix invj(sdim,sdim);
00608                 
00609                 nArray Exi(3), Eeta(3), Ezeta(3);
00610                 clear(Exi); clear(Eeta);
00611                 for ( int i=0; i<nn; ++i ) {
00612                         Exi[0] += coord(i,0)*dNxi(i,0);  // x,xi  
00613                         Exi[1] += coord(i,1)*dNxi(i,0);  // y,xi  
00614                         Exi[2] += coord(i,2)*dNxi(i,0);  // z,xi
00615                         Eeta[0] += coord(i,0)*dNxi(i,1);  // x,eta  
00616                         Eeta[1] += coord(i,1)*dNxi(i,1);  // y,eta 
00617                         Eeta[2] += coord(i,2)*dNxi(i,1);  // z,eta
00618                 }
00619 
00620                 cross_prod(Exi,Eeta,Ezeta);
00621                 invj(0,0)=Exi[0];  invj(0,1)=Eeta[0]; invj(0,2)=Ezeta[0];
00622                 invj(1,0)=Exi[1];  invj(1,1)=Eeta[1]; invj(1,2)=Ezeta[1];
00623                 invj(2,0)=Exi[2];  invj(2,1)=Eeta[2]; invj(2,2)=Ezeta[2];
00624 
00625                 Numeric j=inv3by3(invj);
00626                 for ( int i=0; i<nn; ++i)
00627                 {
00628                         dNx(i,0)=dNxi(i,0)*invj(0,0)+dNxi(i,1)*invj(1,0); // dNIdxi dxidx + dNideta detadx
00629                         dNx(i,1)=dNxi(i,0)*invj(0,1)+dNxi(i,1)*invj(1,1); // dNIdxi dxidy + dNideta detady
00630                         dNx(i,2)=dNxi(i,0)*invj(0,2)+dNxi(i,1)*invj(1,2); // dNIdxi dxidz + dNideta detadz
00631                 }
00632         
00633                 return vect_norm2( Ezeta );
00634         }
00635 
00636         
00643         void cmat_elastic_pstress( Numeric E, Numeric nu, nMatrix &C )
00644         {
00645 #ifdef ALLOW_DYNAMIC_RESIZE
00646                 if ( mat_nrows(C)!=3 || mat_ncols(C)!=3 ) resize( C, 3, 3 );
00647 #elif EXPLICIT_BOUNDS_CHECK 
00648                 if ( mat_nrows(C)!=3 || mat_ncols(C)!=3 ) ErrorMessage("cmat_elastic_pstress","C must be 3 by 3");
00649 #endif
00650                 C(0,0)=E/(1-nu*nu); C(0,1)=nu*C(0,0); C(0,2)=0.0;
00651                 C(1,0)=C(0,1);      C(1,1)=C(0,0);    C(1,2)=0.0;
00652                 C(2,0)=0.0;         C(2,1)=0.0;       C(2,2)=0.5*(1-nu)*C(0,0);
00653         }
00654         
00661         void cmat_elastic_pstrain( Numeric E, Numeric nu, nMatrix &C )
00662         {
00663 #ifdef ALLOW_DYNAMIC_RESIZE
00664                 if ( mat_nrows(C)!=3 || mat_ncols(C)!=3 ) resize( C, 3, 3 );
00665 #elif EXPLICIT_BOUNDS_CHECK 
00666                 if ( mat_nrows(C)!=3 || mat_ncols(C)!=3 ) ErrorMessage("cmat_elastic_pstrain","C must be 3 by 3");
00667 #endif
00668                 Numeric c1, c2, c3;
00669                 c1=E*(1-nu)/(1+nu)/(1-2*nu);
00670                 c2=E*nu/(1+nu)/(1-2*nu);
00671                 c3=0.5*E/(1+nu);
00672                 C(0,0)=c1; C(0,1)=c2; C(0,2)=0.0;  
00673                 C(1,0)=c2; C(1,1)=c1; C(1,2)=0.0;
00674                 C(2,0)=0.0; C(2,1)=0.0; C(2,2)=c3;
00675         }       
00676         
00683         void cmat_elastic_3d( Numeric E, Numeric nu, nMatrix &C )
00684         {
00685 #ifdef ALLOW_DYNAMIC_RESIZE
00686                 if ( mat_nrows(C)!=6 || mat_ncols(C)!=6 ) resize( C, 6, 6 );
00687 #elif EXPLICIT_BOUNDS_CHECK 
00688                 if ( mat_nrows(C)!=6 || mat_ncols(C)!=6 ) ErrorMessage("cmat_elastic_3d","C must be 6 by 6");
00689 #endif
00690                 Numeric c1, c2, c3;
00691                 c1=E*(1-nu)/(1+nu)/(1-2*nu);
00692                 c2=E*nu/(1+nu)/(1-2*nu);
00693                 c3=0.5*E/(1+nu);
00694                 C(0,0)=c1; C(0,1)=c2; C(0,2)=c2; C(0,3)=0.0; C(0,4)=0.0; C(0,5)=0.0;
00695                 C(1,0)=c2; C(1,1)=c1; C(1,2)=c2; C(1,3)=0.0; C(1,4)=0.0; C(1,5)=0.0;
00696                 C(2,0)=c2; C(2,1)=c2; C(2,2)=c1; C(2,3)=0.0; C(2,4)=0.0; C(2,5)=0.0;
00697                 C(3,0)=0.0; C(3,1)=0.0; C(3,2)=0.0; C(3,3)=c3; C(3,4)=0.0; C(3,5)=0.0;
00698                 C(4,0)=0.0; C(4,1)=0.0; C(4,2)=0.0; C(4,3)=0.0; C(4,4)=c3; C(4,5)=0.0;
00699                 C(5,0)=0.0; C(5,1)=0.0; C(5,2)=0.0; C(5,3)=0.0; C(5,4)=0.0; C(5,5)=c3;
00700         }       
00701         
00712         template < class MATA, class MATB, class MATC, class MATK >
00713         void femmult( const MATA &A, const MATB &B, const MATC &C, Numeric alpha, MATK &K )
00714         {
00715                 int n=mat_nrows(A), m=mat_ncols(C), k=mat_ncols(A);
00716 #ifdef ALLOW_DYNAMIC_RESIZE
00717                 if ( n!=mat_nrows(K) || m!=mat_ncols(K) ) resize( K, n, m );
00718 #elif EXPLICIT_BOUNDS_CHECK 
00719                 if ( mat_ncols(A)!=mat_nrows(B) || mat_ncols(B)!=mat_nrows(C) 
00720                         || n!=mat_nrows(K) || m!=mat_ncols(K) ) 
00721                                 ErrorMessage("femmult","Incorrect matrix dimensions");
00722 #endif          
00723                 nMatrix AB(n,k); 
00724                 mult(A,B,AB);
00725                 mult(AB,C,K);
00726                 scale(K,alpha);
00727         }       
00728         
00741         template < class MATA, class MATB, class MATC, class MATK >
00742         void femmult_add( const MATA &A, const MATB &B, const MATC &C, Numeric alpha, MATK &K )
00743         {
00744                 int n=mat_nrows(A), m=mat_ncols(C), k=mat_ncols(A);
00745 #ifdef ALLOW_DYNAMIC_RESIZE
00746                 if ( n!=mat_nrows(K) || m!=mat_ncols(K) ) resize( K, n, m );
00747 #elif EXPLICIT_BOUNDS_CHECK 
00748                 if ( mat_ncols(A)!=mat_nrows(B) || mat_ncols(B)!=mat_nrows(C) 
00749                                 || n!=mat_nrows(K) || m!=mat_ncols(K) ) 
00750                                         ErrorMessage("femmult_add","Incorrect matrix dimensions");
00751 #endif          
00752                 nMatrix AB(n,k); 
00753                 mult(A,B,AB);
00754 //              mult_add(AB,C,K); This does not work for some odd reason it seems to be a gmm bug!
00755                 nMatrix Kadd(n,m);
00756                 mult(AB,C,Kadd);
00757                 add(scaled(Kadd,alpha),K);
00758                 
00759         }
00760 
00772         template < class MATA, class MATB, class MATK >
00773         void femmult_add( const MATA &A, const MATB &B, Numeric alpha, MATK &K )
00774         {
00775                 int n=mat_nrows(A), m=mat_ncols(B);
00776 #ifdef ALLOW_DYNAMIC_RESIZE
00777                 if ( n!=mat_nrows(K) || m!=mat_ncols(K) ) resize( K, n, m );
00778 #elif EXPLICIT_BOUNDS_CHECK 
00779                 if ( mat_ncols(A)!=mat_nrows(B) || n!=mat_nrows(K) || m!=mat_ncols(K) ) 
00780                         ErrorMessage("femmult_add","Incorrect matrix dimensions");
00781 #endif          
00782                 nMatrix AB(n,m); 
00783                 mult(A,B,AB);
00784                 add(scaled(AB,alpha),K);
00785 
00786         }
00787 
00788 
00789                 
00790 } // end namespace
00791 
00792 #endif
00793 
00794
Generated on Sat Jan 19 09:03:57 2008 for Basso by  doxygen 1.5.2