Merge pull request #89 from APSIS-ANALYSIS/luo-dev

Luo dev：使用了名字空间来管理不同张量类的辅助函数。

examples/fsi/src/PLocAssem_2x2Block_Tet4_VMS_Hyperelasticity.cpp

@@ -184,7 +184,7 @@ void PLocAssem_2x2Block_Tet4_VMS_Hyperelasticity::Assem_Residual(
 
     const Tensor2_3D F( ux_x + 1.0, ux_y, ux_z, uy_x, uy_y + 1.0, uy_z, uz_x, uz_y, uz_z + 1.0 );
 
-    const Tensor2_3D invF = inverse( F );
+    const Tensor2_3D invF = Ten2::inverse( F );
 
     const Tensor2_3D DVelo( vx_x, vx_y, vx_z, vy_x, vy_y, vy_z, vz_x, vz_y, vz_z );
 
@@ -200,7 +200,7 @@ void PLocAssem_2x2Block_Tet4_VMS_Hyperelasticity::Assem_Residual(
         qua_prestress[qua*6+5], qua_prestress[qua*6+1], qua_prestress[qua*6+3],
         qua_prestress[qua*6+4], qua_prestress[qua*6+3], qua_prestress[qua*6+2] );
 
-    P_iso += prestress * cofactor( F );
+    P_iso += prestress * Ten2::cofactor( F );
     // ------------------------------------------------------------------------
 
     const double rho = matmodel->get_rho(p);
@@ -343,7 +343,7 @@ void PLocAssem_2x2Block_Tet4_VMS_Hyperelasticity::Assem_Tangent_Residual(
 
     const Tensor2_3D F( ux_x + 1.0, ux_y, ux_z, uy_x, uy_y + 1.0, uy_z, uz_x, uz_y, uz_z + 1.0 );
 
-    const Tensor2_3D invF = inverse( F );
+    const Tensor2_3D invF = Ten2::inverse( F );
 
     const Tensor2_3D DVelo( vx_x, vx_y, vx_z, vy_x, vy_y, vy_z, vz_x, vz_y, vz_z );
 
@@ -364,7 +364,7 @@ void PLocAssem_2x2Block_Tet4_VMS_Hyperelasticity::Assem_Tangent_Residual(
         qua_prestress[qua*6+5], qua_prestress[qua*6+1], qua_prestress[qua*6+3],
         qua_prestress[qua*6+4], qua_prestress[qua*6+3], qua_prestress[qua*6+2] );
 
-    P_iso += prestress * cofactor( F );
+    P_iso += prestress * Ten2::cofactor( F );
     // ------------------------------------------------------------------------
 
     const double rho = matmodel->get_rho(p);
@@ -635,7 +635,7 @@ void PLocAssem_2x2Block_Tet4_VMS_Hyperelasticity::Assem_Mass_Residual(
 
     const Tensor2_3D F( ux_x + 1.0, ux_y, ux_z, uy_x, uy_y + 1.0, uy_z, uz_x, uz_y, uz_z + 1.0 );
 
-    const Tensor2_3D invF = inverse( F );
+    const Tensor2_3D invF = Ten2::inverse( F );
 
     const Tensor2_3D DVelo( vx_x, vx_y, vx_z, vy_x, vy_y, vy_z, vz_x, vz_y, vz_z );
 
@@ -651,7 +651,7 @@ void PLocAssem_2x2Block_Tet4_VMS_Hyperelasticity::Assem_Mass_Residual(
         qua_prestress[qua*6+5], qua_prestress[qua*6+1], qua_prestress[qua*6+3],
         qua_prestress[qua*6+4], qua_prestress[qua*6+3], qua_prestress[qua*6+2] );
 
-    P_iso += prestress * cofactor( F );
+    P_iso += prestress * Ten2::cofactor( F );
     // ------------------------------------------------------------------------
 
     double mbeta = matmodel->get_beta(p);

examples/fsi/src/PLocAssem_2x2Block_Tet4_VMS_Incompressible.cpp

@@ -192,7 +192,7 @@ void PLocAssem_2x2Block_Tet4_VMS_Incompressible::Assem_Residual(
 
     const Tensor2_3D F( ux_x + 1.0, ux_y, ux_z, uy_x, uy_y + 1.0, uy_z, uz_x, uz_y, uz_z + 1.0 );
 
-    const Tensor2_3D invF = inverse(F);
+    const Tensor2_3D invF = Ten2::inverse(F);
 
     const Tensor2_3D DVelo( vx_x, vx_y, vx_z, vy_x, vy_y, vy_z, vz_x, vz_y, vz_z );
 
@@ -207,7 +207,7 @@ void PLocAssem_2x2Block_Tet4_VMS_Incompressible::Assem_Residual(
         qua_prestress[qua*6+5], qua_prestress[qua*6+1], qua_prestress[qua*6+3],
         qua_prestress[qua*6+4], qua_prestress[qua*6+3], qua_prestress[qua*6+2] );
 
-    P_iso += prestress * cofactor( F );
+    P_iso += prestress * Ten2::cofactor( F );
     // ------------------------------------------------------------------------
 
     const double rho = matmodel->get_rho(p);
@@ -359,7 +359,7 @@ void PLocAssem_2x2Block_Tet4_VMS_Incompressible::Assem_Tangent_Residual(
 
     const Tensor2_3D F( ux_x + 1.0, ux_y, ux_z, uy_x, uy_y + 1.0, uy_z, uz_x, uz_y, uz_z + 1.0 );
 
-    const Tensor2_3D invF = inverse(F);
+    const Tensor2_3D invF = Ten2::inverse(F);
 
     const Tensor2_3D DVelo( vx_x, vx_y, vx_z, vy_x, vy_y, vy_z, vz_x, vz_y, vz_z );
 
@@ -380,7 +380,7 @@ void PLocAssem_2x2Block_Tet4_VMS_Incompressible::Assem_Tangent_Residual(
         qua_prestress[qua*6+5], qua_prestress[qua*6+1], qua_prestress[qua*6+3],
         qua_prestress[qua*6+4], qua_prestress[qua*6+3], qua_prestress[qua*6+2] );
 
-    P_iso += prestress * cofactor( F );
+    P_iso += prestress * Ten2::cofactor( F );
     // ------------------------------------------------------------------------
 
     const double rho = matmodel->get_rho(p);
@@ -629,7 +629,7 @@ void PLocAssem_2x2Block_Tet4_VMS_Incompressible::Assem_Mass_Residual(
 
     const Tensor2_3D F( ux_x + 1.0, ux_y, ux_z, uy_x, uy_y + 1.0, uy_z, uz_x, uz_y, uz_z + 1.0 );
 
-    const Tensor2_3D invF = inverse(F);
+    const Tensor2_3D invF = Ten2::inverse(F);
 
     const Tensor2_3D DVelo(  vx_x, vx_y, vx_z, vy_x, vy_y, vy_z, vz_x, vz_y, vz_z );
 
@@ -645,7 +645,7 @@ void PLocAssem_2x2Block_Tet4_VMS_Incompressible::Assem_Mass_Residual(
         qua_prestress[qua*6+5], qua_prestress[qua*6+1], qua_prestress[qua*6+3],
         qua_prestress[qua*6+4], qua_prestress[qua*6+3], qua_prestress[qua*6+2] );
 
-    P_iso += prestress * cofactor( F );
+    P_iso += prestress * Ten2::cofactor( F );
     // ------------------------------------------------------------------------
 
     double mbeta = matmodel->get_beta(p);

examples/fsi/vis_fsi_wss.cpp

@@ -181,7 +181,7 @@ int main( int argc, char * argv[] )
         v_ctrlPts[3*trn[2]+1] - v_ctrlPts[3*trn[0]+1],
         v_ctrlPts[3*trn[2]+2] - v_ctrlPts[3*trn[0]+2] );
 
-    outnormal[ee] = VEC3_T::cross_product(vec01, vec02); // out = vec01 x vec02
+    outnormal[ee] = Vec3::cross_product(vec01, vec02); // out = vec01 x vec02
 
     // return out length and scale itself to have unit length
     tri_area[ee] = 0.5 * outnormal[ee].normalize();
@@ -190,7 +190,7 @@ int main( int argc, char * argv[] )
         v_ctrlPts[interior_node[ee]*3+1] - v_ctrlPts[3*trn[0]+1],
         v_ctrlPts[interior_node[ee]*3+2] - v_ctrlPts[3*trn[0]+2] );
 
-    if( VEC3_T::dot_product(outnormal[ee], vec03)> 0 ) outnormal[ee] *= -1.0;
+    if( Vec3::dot_product(outnormal[ee], vec03)> 0 ) outnormal[ee] *= -1.0;
   }
 
   // Clean the volumetric data to save memory

examples/ns/vis_p1_wss.cpp

@@ -165,15 +165,15 @@ int main( int argc, char * argv[] )
         v_ctrlPts[3*trn[2]+1] - v_ctrlPts[3*trn[0]+1],
         v_ctrlPts[3*trn[2]+2] - v_ctrlPts[3*trn[0]+2] ); 
 
-    Vector_3 ou = VEC3_T::cross_product( l01, l02 );
+    Vector_3 ou = Vec3::cross_product( l01, l02 );
 
     tri_area[ee] = 0.5 * ou.normalize();
 
     const Vector_3 inw( v_ctrlPts[interior_node[ee]*3]   - v_ctrlPts[3*trn[0]],
         v_ctrlPts[interior_node[ee]*3+1] - v_ctrlPts[3*trn[0]+1],
         v_ctrlPts[interior_node[ee]*3+2] - v_ctrlPts[3*trn[0]+2] );
 
-    const double out_dot_in = VEC3_T::dot_product( ou, inw );
+    const double out_dot_in = Vec3::dot_product( ou, inw );
 
     // if in case the normal points inside, multiply by -1
     if( out_dot_in > 0 ) ou *= -1.0;

examples/ruc_fsi/src/PGAssem_Tet_Wall.cpp

@@ -187,7 +187,7 @@ void PGAssem_Tet_Wall::RingBC_KG(
         const Tensor2_3D QT = ringnbc_part->get_rotation_matrix( pos );
 
         // Skew-to-global transformation matrix
-        const Tensor2_3D Q = transpose( QT );
+        const Tensor2_3D Q = Ten2::transpose( QT );
 
         for( int jj = dof-1; jj < ncol; jj += dof )
         {

examples/ruc_fsi/vis_p1_wss.cpp

@@ -167,15 +167,15 @@ int main( int argc, char * argv[] )
      v_ctrlPts[3*trn[2]+1] - v_ctrlPts[3*trn[0]+1],
      v_ctrlPts[3*trn[2]+2] - v_ctrlPts[3*trn[0]+2] ); 
 
-    Vector_3 ou = VEC3_T::cross_product( l01, l02 );
+    Vector_3 ou = Vec3::cross_product( l01, l02 );
 
     tri_area[ee] = 0.5 * ou.normalize();
 
     const Vector_3 inw( v_ctrlPts[interior_node[ee]*3]   - v_ctrlPts[3*trn[0]],
         v_ctrlPts[interior_node[ee]*3+1] - v_ctrlPts[3*trn[0]+1],
         v_ctrlPts[interior_node[ee]*3+2] - v_ctrlPts[3*trn[0]+2] );
 
-    const double out_dot_in = VEC3_T::dot_product( ou, inw ); 
+    const double out_dot_in = Vec3::dot_product( ou, inw ); 
 
     // if in case the normal points inside, multiply by -1
     if(out_dot_in > 0) ou *= -1.0;

examples/test/elem_test.cpp

@@ -160,15 +160,15 @@ int main( int argc, char * argv[] )
   vec1.gen_rand();
   vec2.gen_rand();
   vec3.gen_rand();
-  Vector_3 vec_cross = VEC3_T::cross_product(vec1, vec2);
+  Vector_3 vec_cross = Vec3::cross_product(vec1, vec2);
   if(vec_cross.dot_product( vec3 )<0.0)
   {
     Vector_3 temp = vec1;
     vec1 = vec2;
     vec2 = temp;
   }
 
-  const double volume = vec3.dot_product(VEC3_T::cross_product(vec1, vec2));
+  const double volume = vec3.dot_product(Vec3::cross_product(vec1, vec2));
 
   double * ctrl_x_hex = new double[8]{0.0, vec1.x(), vec1.x()+vec2.x(), vec2.x(), 
     vec3.x(), vec3.x()+vec1.x(), vec3.x()+vec1.x()+vec2.x(), vec3.x()+vec2.x()};

include/IMaterialModel.hpp

@@ -93,7 +93,7 @@ class IMaterialModel
     {
       Tensor2_3D P, S;
       get_PK(F, P, S);
-      const Tensor2_3D Ft = transpose( F );
+      const Tensor2_3D Ft = Ten2::transpose( F );
       Tensor2_3D sigma = P * Ft;
       sigma.scale( (1.0/F.det()) );
       return sigma;

include/SymmTensor2_3D.hpp

@@ -197,17 +197,20 @@ Tensor2_3D operator*( const SymmTensor2_3D &left, const SymmTensor2_3D &right );
 
 SymmTensor2_3D operator*( const double &val, const SymmTensor2_3D &input );
 
-// Return the inverse of the input matrix
-SymmTensor2_3D inverse( const SymmTensor2_3D &input );
+namespace STen2
+{
+  // Return the inverse of the input matrix
+  SymmTensor2_3D inverse( const SymmTensor2_3D &input );
 
-// Generate the right Cauchy-Green tensor C = F^T F
-SymmTensor2_3D gen_right_Cauchy_Green( const Tensor2_3D &input );
+  // Generate the right Cauchy-Green tensor C = F^T F
+  SymmTensor2_3D gen_right_Cauchy_Green( const Tensor2_3D &input );
 
-// Generate the left Cauchy-Green tensor b = F F^T
-SymmTensor2_3D gen_left_Cauchy_Green( const Tensor2_3D &input );
+  // Generate the left Cauchy-Green tensor b = F F^T
+  SymmTensor2_3D gen_left_Cauchy_Green( const Tensor2_3D &input );
 
-// Convert a regular matrix to its symmetric part
-// output = 0.5 x ( source + source_transpose )
-SymmTensor2_3D gen_symm_part( const Tensor2_3D &input );
+  // Convert a regular matrix to its symmetric part
+  // output = 0.5 x ( source + source_transpose )
+  SymmTensor2_3D gen_symm_part( const Tensor2_3D &input );
+}
 
 #endif

include/SymmTensor4_3D.hpp

@@ -217,12 +217,14 @@ class SymmTensor4_3D
 };
 
 // These functions behave in an identical manner to the member function of the
-// same function name. For the gen_zero and gem_symm_id, we added a middle name
-// ST4 to differentiate them from the functions of Tensor4_3D.
-SymmTensor4_3D gen_ST4_zero();
+// same function name.
+namespace STen4
+{
+  SymmTensor4_3D gen_zero();
 
-SymmTensor4_3D gen_ST4_symm_id();
+  SymmTensor4_3D gen_symm_id();
 
-SymmTensor4_3D gen_ST4_Ptilde( const SymmTensor2_3D &invC );
+  SymmTensor4_3D gen_Ptilde( const SymmTensor2_3D &invC );
+}
 
 #endif

include/Tensor2_3D.hpp

@@ -270,19 +270,22 @@ Tensor2_3D operator*( const Tensor2_3D &left, const Tensor2_3D &right );
 // Return the scalar scaling of matrix
 Tensor2_3D operator*( const double &val, const Tensor2_3D &input );
 
-// Return the inverse of the input matrix
-Tensor2_3D inverse( const Tensor2_3D &input );
+namespace Ten2
+{
+  // Return the inverse of the input matrix
+  Tensor2_3D inverse( const Tensor2_3D &input );
 
-// Return the cofactor of input matrix which is J input^(-T)
-Tensor2_3D cofactor( const Tensor2_3D &input );
+  // Return the cofactor of input matrix which is J input^(-T)
+  Tensor2_3D cofactor( const Tensor2_3D &input );
 
-// Return the transpose of input matrix
-Tensor2_3D transpose( const Tensor2_3D &input );
+  // Return the transpose of input matrix
+  Tensor2_3D transpose( const Tensor2_3D &input );
 
-// Return an identity matrix
-Tensor2_3D gen_identity_matrix();
+  // Return an identity matrix
+  Tensor2_3D gen_id();
 
-// Return a zero matrix
-Tensor2_3D gen_zero_matrix();
+  // Return a zero matrix
+  Tensor2_3D gen_zero();
+}
 
 #endif

include/Tensor4_3D.hpp

@@ -128,7 +128,7 @@ class Tensor4_3D
     // ------------------------------------------------------------------------
     void gen_P( const Tensor2_3D &C, const Tensor2_3D &invC );
 
-    void gen_P( const Tensor2_3D &C ) { gen_P( C, inverse(C) ); }
+    void gen_P( const Tensor2_3D &C ) { gen_P( C, Ten2::inverse(C) ); }
 
     // ------------------------------------------------------------------------
     // Generate Projector Ptilde = invC O invC - 1/3 invC x invC
@@ -302,35 +302,38 @@ Tensor4_3D operator*( const Tensor4_3D &tleft, const Tensor4_3D &tright );
 // Return scalar multiplication on the input tensor
 Tensor4_3D operator*( const double &val, const Tensor4_3D &input );
 
-Tensor4_3D gen_T4_zero();
-
-Tensor4_3D gen_T4_symm_id();
-
-// ------------------------------------------------------------------------
-// Generate Projector P = SymmId4 - 1/3 invC x C
-// P_IJKL = SymmID_IJKL - 1/3 invC_IJ C_KL
-// C is assumed to be the right Cauchy-Green tensor
-// invC is the inverse of C
-// see Holzapfel book p.229 eqn. (6.84).
-// ------------------------------------------------------------------------
-Tensor4_3D gen_T4_P( const Tensor2_3D &C, const Tensor2_3D &invC );
-
-Tensor4_3D gen_T4_P( const Tensor2_3D &C );
-
-// ------------------------------------------------------------------------
-// Generate Projector Pt = transpose of P = SymmId4 - 1/3 C x invC
-// P_IJKL = SymmID_IJKL - 1/3 C_IJ invC_KL
-// C is assumed to be the right Cauchy-Green tensor
-// invC is the inverse of C
-// see Holzapfel book p.229 eqn. (6.84).
-// ------------------------------------------------------------------------
-Tensor4_3D gen_T4_Pt( const Tensor2_3D &C );
-
-// ------------------------------------------------------------------------
-// Generate Projector Ptilde = invC O invC - 1/3 invC x invC
-// invC is assumed to be the right Cauchy-Green tensor 
-// see Holzapfel book p. 255, eqn. (6.170).
-// ------------------------------------------------------------------------
-Tensor4_3D gen_T4_Ptilde( const Tensor2_3D &invC );
+namespace Ten4
+{
+  Tensor4_3D gen_zero();
+
+  Tensor4_3D gen_symm_id();
+
+  // ------------------------------------------------------------------------
+  // Generate Projector P = SymmId4 - 1/3 invC x C
+  // P_IJKL = SymmID_IJKL - 1/3 invC_IJ C_KL
+  // C is assumed to be the right Cauchy-Green tensor
+  // invC is the inverse of C
+  // see Holzapfel book p.229 eqn. (6.84).
+  // ------------------------------------------------------------------------
+  Tensor4_3D gen_P( const Tensor2_3D &C, const Tensor2_3D &invC );
+
+  Tensor4_3D gen_P( const Tensor2_3D &C );
+
+  // ------------------------------------------------------------------------
+  // Generate Projector Pt = transpose of P = SymmId4 - 1/3 C x invC
+  // P_IJKL = SymmID_IJKL - 1/3 C_IJ invC_KL
+  // C is assumed to be the right Cauchy-Green tensor
+  // invC is the inverse of C
+  // see Holzapfel book p.229 eqn. (6.84).
+  // ------------------------------------------------------------------------
+  Tensor4_3D gen_Pt( const Tensor2_3D &C );
+
+  // ------------------------------------------------------------------------
+  // Generate Projector Ptilde = invC O invC - 1/3 invC x invC
+  // invC is assumed to be the right Cauchy-Green tensor 
+  // see Holzapfel book p. 255, eqn. (6.170).
+  // ------------------------------------------------------------------------
+  Tensor4_3D gen_Ptilde( const Tensor2_3D &invC );
+}
 
 #endif

include/Vector_3.hpp

@@ -105,7 +105,7 @@ class Vector_3
 // calculate a scalar product of a input vector
 Vector_3 operator*( const double &val, const Vector_3 &source );
 
-namespace VEC3_T
+namespace Vec3
 {
   // calculate the distance between two vector by L2 norm
   double dist( const Vector_3 &a, const Vector_3 &b );

src/Analysis_Tool/ALocal_InflowBC.cpp

@@ -111,14 +111,14 @@ double ALocal_InflowBC::get_radius( const int &nbc_id,
   // inflow boundary points (i.e. Num_LD = 0 ).
   SYS_T::print_fatal_if( num_out_bc_pts[nbc_id] == 0, "Error: ALocal_InflowBC::get_radius, this function can only be called in sub-domains which contains the inflow boundary node.\n");
 
-  const double rc = VEC3_T::dist( pt, centroid[nbc_id] );
+  const double rc = Vec3::dist( pt, centroid[nbc_id] );
 
   // Now loop over the boundary points to find rb.
-  double rb = VEC3_T::dist( pt, Vector_3( outline_pts[nbc_id][0], outline_pts[nbc_id][1], outline_pts[nbc_id][2]) );
+  double rb = Vec3::dist( pt, Vector_3( outline_pts[nbc_id][0], outline_pts[nbc_id][1], outline_pts[nbc_id][2]) );
 
   for(int ii=1; ii<num_out_bc_pts[nbc_id]; ++ii)
   {
-    double newdist = VEC3_T::dist( pt, Vector_3( outline_pts[nbc_id][3*ii], outline_pts[nbc_id][3*ii+1], outline_pts[nbc_id][3*ii+2]) );
+    double newdist = Vec3::dist( pt, Vector_3( outline_pts[nbc_id][3*ii], outline_pts[nbc_id][3*ii+1], outline_pts[nbc_id][3*ii+2]) );
 
     if(newdist < rb) rb = newdist;
   }