trdotn test

neams-th-coe · Dec 11, 2023 · 9109e82 · 9109e82
1 parent e07e0d7
commit 9109e82
Show file tree

Hide file tree

Showing 4 changed files with 36 additions and 22 deletions.
diff --git a/include/base/CardinalEnums.h b/include/base/CardinalEnums.h
@@ -98,7 +98,7 @@ enum NekFieldEnum
   scalar02,
   scalar03,
   wall_shear,
-//  traction,
+  traction,
   traction_x,
   traction_y,
   traction_z,

diff --git a/src/base/CardinalEnums.C b/src/base/CardinalEnums.C
@@ -47,7 +47,7 @@ getNekFieldEnum()
 {
   return MooseEnum(
       "velocity_component velocity_x velocity_y velocity_z velocity temperature"
-      " pressure scalar01 scalar02 scalar03 wall_shear traction_x traction_y traction_z"
+      " pressure scalar01 scalar02 scalar03 wall_shear traction traction_x traction_y traction_z"
       " ros_s11 ros_s22 ros_s33 ros_s12 ros_s23 ros_s13 unity");
 }
 

diff --git a/src/base/NekInterface.C b/src/base/NekInterface.C
@@ -665,7 +665,7 @@ initializeRateOfStrainTensor()
 void
 initializeTraction()
 {
-  nekrs::_traction = (double *)calloc(3*velocityFieldOffset(), sizeof(double));
+  nekrs::_traction = (double *)calloc(4*velocityFieldOffset(), sizeof(double)); //TODO: finalize size
 }
 
 void
@@ -1509,18 +1509,32 @@ computeTraction(double * traction, const nek_mesh::NekMeshEnum pp_mesh)
           int vol_id = mesh->vmapM[offset + v];
           int surf_offset = mesh->Nsgeo * (offset + v);
 
+          double nx =  sgeo[surf_offset + NXID];
+          double ny =  sgeo[surf_offset + NYID];
+          double nz =  sgeo[surf_offset + NZID];
           //TODO: explain reasoning below
-          traction[nrs_offset*0 + vol_id] = Tau_ij[0*nrs_offset + vol_id] * sgeo[surf_offset + NXID]
-                                          + Tau_ij[3*nrs_offset + vol_id] * sgeo[surf_offset + NYID]
-                                          + Tau_ij[5*nrs_offset + vol_id] * sgeo[surf_offset + NZID];
+          traction[nrs_offset*0 + vol_id] = Tau_ij[0*nrs_offset + vol_id] * nx
+                                          + Tau_ij[3*nrs_offset + vol_id] * ny
+                                          + Tau_ij[5*nrs_offset + vol_id] * nz;
 
-          traction[nrs_offset*1 + vol_id] = Tau_ij[3*nrs_offset + vol_id] * sgeo[surf_offset + NXID]
-                                          + Tau_ij[1*nrs_offset + vol_id] * sgeo[surf_offset + NYID]
-                                          + Tau_ij[4*nrs_offset + vol_id] * sgeo[surf_offset + NZID];
+          traction[nrs_offset*1 + vol_id] = Tau_ij[3*nrs_offset + vol_id] * nx
+                                          + Tau_ij[1*nrs_offset + vol_id] * ny
+                                          + Tau_ij[4*nrs_offset + vol_id] * nz;
 
-          traction[nrs_offset*2 + vol_id] = Tau_ij[5*nrs_offset + vol_id] * sgeo[surf_offset + NXID]
-                                          + Tau_ij[4*nrs_offset + vol_id] * sgeo[surf_offset + NYID]
-                                          + Tau_ij[2*nrs_offset + vol_id] * sgeo[surf_offset + NZID];
+          traction[nrs_offset*2 + vol_id] = Tau_ij[5*nrs_offset + vol_id] * nx
+                                          + Tau_ij[4*nrs_offset + vol_id] * ny
+                                          + Tau_ij[2*nrs_offset + vol_id] * nz;
+
+          double t_dot_n = traction[0*nrs_offset + vol_id] * nx  // TODO: decide whether to keep
+                         + traction[1*nrs_offset + vol_id] * ny
+                         + traction[2*nrs_offset + vol_id] * nz;
+
+          double tx = traction[nrs_offset*0 + vol_id];
+          double ty = traction[nrs_offset*1 + vol_id];
+          double tz = traction[nrs_offset*2 + vol_id];
+
+          traction[nrs_offset*3 + vol_id] = (t_dot_n/t_dot_n) *
+                                          std::sqrt( tx * tx + ty * ty + tz * tz);
         }
       }
     }
@@ -1758,7 +1772,7 @@ double
 traction(const int id)
 {
   nrs_t * nrs = (nrs_t *)nrsPtr();
-  return _traction[id];
+  return _traction[id+ 3 * nrs->fieldOffset]; //TODO: decide whether to keep this or not
 }
 
 double
@@ -1933,9 +1947,9 @@ double (*solutionPointer(const field::NekFieldEnum & field))(int)
     case field::wall_shear:
       f = &wall_shear;
       break;
-//    case field::traction:
-//      f = &traction;
-//      break;
+    case field::traction:
+      f = &traction;
+      break;
     case field::traction_x:
       f = &traction_x;
       break;
@@ -2102,9 +2116,9 @@ dimensionalize(const field::NekFieldEnum & field, double & value)
     case field::wall_shear:
       // TODO: add dimensionalization
       break;
-//    case field::traction:
-//      // TODO: add dimensionalization
-//      break;
+    case field::traction:
+      // TODO: add dimensionalization
+      break;
     case field::traction_x:
       // TODO: add dimensionalization
       break;

diff --git a/src/base/NekRSProblemBase.C b/src/base/NekRSProblemBase.C
@@ -912,8 +912,8 @@ NekRSProblemBase::extractOutputs()
         field_enum = field::scalar03;
       else if (_var_names[i] == "wall_shear")
         field_enum = field::wall_shear;
-//      else if (_var_names[i] == "traction")
-//        field_enum = field::traction;
+      else if (_var_names[i] == "traction")
+        field_enum = field::traction;
       else if (_var_names[i] == "traction_x")
         field_enum = field::traction_x;
       else if (_var_names[i] == "traction_y")
@@ -1013,7 +1013,7 @@ NekRSProblemBase::addExternalVariables()
       }
       else if (output == "traction")
       {
-//        _var_names.push_back("traction");
+        _var_names.push_back("traction");
         _var_names.push_back("traction_x");
         _var_names.push_back("traction_y");
         _var_names.push_back("traction_z");