diff --git a/src/constitutive/TACSBasicBeamConstitutive.cpp b/src/constitutive/TACSBasicBeamConstitutive.cpp
index 61b525d59..3f740e727 100644
--- a/src/constitutive/TACSBasicBeamConstitutive.cpp
+++ b/src/constitutive/TACSBasicBeamConstitutive.cpp
@@ -122,7 +122,7 @@ TACSBasicBeamConstitutive::TACSBasicBeamConstitutive(
   C[0] = E * A;
   C[6] = G * J;
   C[11] = E * Iz;
-  C[12] = E * Iyz;
+  C[12] = -E * Iyz;
   C[15] = E * Iy;
   C[18] = ky * G * A;
   C[20] = kz * G * A;
@@ -131,9 +131,9 @@ TACSBasicBeamConstitutive::TACSBasicBeamConstitutive(
   rho[0] = density * A;
   rho[1] = 0.0;
   rho[2] = 0.0;
-  rho[3] = density * Iy;
-  rho[4] = density * Iz;
-  rho[5] = density * Iyz;
+  rho[3] = density * Iz;
+  rho[4] = density * Iy;
+  rho[5] = -density * Iyz;
 }
 
 TACSBasicBeamConstitutive::~TACSBasicBeamConstitutive() {
diff --git a/src/constitutive/TACSCompositeShellConstitutive.cpp b/src/constitutive/TACSCompositeShellConstitutive.cpp
index 5082987e4..06606b624 100644
--- a/src/constitutive/TACSCompositeShellConstitutive.cpp
+++ b/src/constitutive/TACSCompositeShellConstitutive.cpp
@@ -27,7 +27,7 @@ const char *TACSCompositeShellConstitutive::constName =
 TACSCompositeShellConstitutive::TACSCompositeShellConstitutive(
     int _num_plies, TACSOrthotropicPly **_ply_props,
     const TacsScalar *_ply_thickness, const TacsScalar *_ply_angles,
-    TacsScalar _kcorr) {
+    TacsScalar _kcorr, TacsScalar _tOffset) {
   num_plies = _num_plies;
   ply_thickness = new TacsScalar[num_plies];
   ply_angles = new TacsScalar[num_plies];
@@ -42,6 +42,7 @@ TACSCompositeShellConstitutive::TACSCompositeShellConstitutive(
   }
 
   kcorr = _kcorr;
+  tOffset = _tOffset;
 }
 
 TACSCompositeShellConstitutive::~TACSCompositeShellConstitutive() {
@@ -83,7 +84,7 @@ void TACSCompositeShellConstitutive::evalMassMoments(int elemIndex,
   }
 
   // Compute the contribution to the mass moment from each layer
-  TacsScalar t0 = -0.5 * t;
+  TacsScalar t0 = -(0.5 + tOffset) * t;
   for (int i = 0; i < num_plies; i++) {
     TacsScalar rho_ply = ply_props[i]->getDensity();
     TacsScalar t1 = t0 + ply_thickness[i];
@@ -131,7 +132,7 @@ void TACSCompositeShellConstitutive::evalStress(int elemIndex,
   }
 
   // Compute the contribution to the stiffness from each layer
-  TacsScalar t0 = -0.5 * t;
+  TacsScalar t0 = -(0.5 + tOffset) * t;
   for (int k = 0; k < num_plies; k++) {
     TacsScalar Qbar[6], Abar[3];
     ply_props[k]->calculateQbar(ply_angles[k], Qbar);
@@ -174,7 +175,7 @@ TacsScalar TACSCompositeShellConstitutive::evalFailure(
   for (int i = 0; i < num_plies; i++) {
     t += ply_thickness[i];
   }
-  TacsScalar t0 = -0.5 * t;
+  TacsScalar t0 = -(0.5 + tOffset) * t;
 
   // Keep track of the maximum failure criterion
   TacsScalar max = 0.0;
@@ -207,7 +208,7 @@ TacsScalar TACSCompositeShellConstitutive::evalFailureStrainSens(
   for (int i = 0; i < num_plies; i++) {
     t += ply_thickness[i];
   }
-  TacsScalar t0 = -0.5 * t;
+  TacsScalar t0 = -(0.5 + tOffset) * t;
 
   // Keep track of the maximum failure criterion
   TacsScalar max = 0.0;
@@ -270,7 +271,7 @@ void TACSCompositeShellConstitutive::evalTangentStiffness(int elemIndex,
   }
 
   // Compute the contribution to the stiffness from each layer
-  TacsScalar t0 = -0.5 * t;
+  TacsScalar t0 = -(0.5 + tOffset) * t;
   for (int k = 0; k < num_plies; k++) {
     TacsScalar Qbar[6], Abar[3];
     ply_props[k]->calculateQbar(ply_angles[k], Qbar);
@@ -353,3 +354,10 @@ void TACSCompositeShellConstitutive::getPlyAngles(TacsScalar *_ply_angles) {
     _ply_angles[i] = ply_angles[i];
   }
 }
+
+/*
+  Get thickness offset
+*/
+TacsScalar TACSCompositeShellConstitutive::getThicknessOffset() {
+  return tOffset;
+}
diff --git a/src/constitutive/TACSCompositeShellConstitutive.h b/src/constitutive/TACSCompositeShellConstitutive.h
index 42cd0068a..268354294 100644
--- a/src/constitutive/TACSCompositeShellConstitutive.h
+++ b/src/constitutive/TACSCompositeShellConstitutive.h
@@ -28,7 +28,8 @@ class TACSCompositeShellConstitutive : public TACSShellConstitutive {
                                  TACSOrthotropicPly **_ply_props,
                                  const TacsScalar *_ply_thickness,
                                  const TacsScalar *_ply_angles,
-                                 TacsScalar _kcorr = 5.0 / 6.0);
+                                 TacsScalar _kcorr = 5.0 / 6.0,
+                                 TacsScalar _tOffset = 0.0);
   ~TACSCompositeShellConstitutive();
 
   // Evaluate the material density
@@ -78,13 +79,14 @@ class TACSCompositeShellConstitutive : public TACSShellConstitutive {
   // Get ply angles and thicknesses
   void getPlyThicknesses(TacsScalar *_ply_thickness);
   void getPlyAngles(TacsScalar *_ply_angles);
+  TacsScalar getThicknessOffset();
 
  private:
   // Store information about the design variable
   int num_plies;
   TACSOrthotropicPly **ply_props;
   TacsScalar *ply_thickness, *ply_angles;
-  TacsScalar kcorr;
+  TacsScalar kcorr, tOffset;
 
   // The object name
   static const char *constName;
diff --git a/src/constitutive/TACSIsoRectangleBeamConstitutive.cpp b/src/constitutive/TACSIsoRectangleBeamConstitutive.cpp
index 950ae9d0c..9f6df6f1d 100644
--- a/src/constitutive/TACSIsoRectangleBeamConstitutive.cpp
+++ b/src/constitutive/TACSIsoRectangleBeamConstitutive.cpp
@@ -4,7 +4,7 @@ TACSIsoRectangleBeamConstitutive::TACSIsoRectangleBeamConstitutive(
     TACSMaterialProperties* properties, TacsScalar _width,
     TacsScalar _thickness, int _width_num, int _thickness_num,
     TacsScalar _lb_width, TacsScalar _ub_width, TacsScalar _lb_thickness,
-    TacsScalar _ub_thickness) {
+    TacsScalar _ub_thickness, TacsScalar _w_offset, TacsScalar _t_offset) {
   props = properties;
   props->incref();
 
@@ -16,6 +16,8 @@ TACSIsoRectangleBeamConstitutive::TACSIsoRectangleBeamConstitutive(
   ub_thickness = _ub_thickness;
   lb_width = _lb_width;
   ub_width = _ub_width;
+  t_offset = _t_offset;
+  w_offset = _w_offset;
 
   ks_weight = 100.0;
 }
@@ -94,17 +96,20 @@ void TACSIsoRectangleBeamConstitutive::evalMassMoments(int elemIndex,
                                                        TacsScalar moments[]) {
   TacsScalar rho = props->getDensity();
   TacsScalar A = thickness * width;
+  TacsScalar delta_y = t_offset * thickness;
+  TacsScalar delta_z = w_offset * width;
   TacsScalar t3 = thickness * thickness * thickness;
   TacsScalar w3 = width * width * width;
-  TacsScalar Iy = 1.0 / 12.0 * thickness * w3;
-  TacsScalar Iz = 1.0 / 12.0 * width * t3;
+  TacsScalar Iy = 1.0 / 12.0 * thickness * w3 + delta_z * delta_z * A;
+  TacsScalar Iz = 1.0 / 12.0 * width * t3 + delta_y * delta_y * A;
+  TacsScalar Iyz = -delta_y * delta_z * A;
 
   moments[0] = rho * A;
-  moments[1] = 0.0;  // centroid offset y?
-  moments[2] = 0.0;  // centroid offset z?
-  moments[3] = rho * Iy;
-  moments[4] = rho * Iz;
-  moments[5] = 0.0;  // rho * Iyz
+  moments[1] = rho * delta_y * A;  // centroid offset y?
+  moments[2] = rho * delta_z * A;  // centroid offset z?
+  moments[3] = rho * Iz;
+  moments[4] = rho * Iy;
+  moments[5] = -rho * Iyz;  // -rho * Iyz
 }
 
 void TACSIsoRectangleBeamConstitutive::addMassMomentsDVSens(
@@ -112,26 +117,39 @@ void TACSIsoRectangleBeamConstitutive::addMassMomentsDVSens(
     const TacsScalar scale[], int dvLen, TacsScalar dfdx[]) {
   TacsScalar rho = props->getDensity();
   TacsScalar A = thickness * width;
+  TacsScalar delta_y = t_offset * thickness;
+  TacsScalar delta_z = w_offset * width;
   TacsScalar t3 = thickness * thickness * thickness;
   TacsScalar w3 = width * width * width;
-  TacsScalar Iy = 1.0 / 12.0 * thickness * w3;
-  TacsScalar Iz = 1.0 / 12.0 * width * t3;
+  TacsScalar Iy = 1.0 / 12.0 * thickness * w3 + delta_z * delta_z * A;
+  TacsScalar Iz = 1.0 / 12.0 * width * t3 + delta_y * delta_y * A;
+  TacsScalar Iyz = -delta_y * delta_z * A;
 
   int index = 0;
   if (width_num >= 0) {
     TacsScalar dA = A / width;
-    TacsScalar dIy = 3.0 * Iy / width;
-    TacsScalar dIz = Iz / width;
-
-    dfdx[index] += rho * (scale[0] * dA + scale[3] * dIy + scale[4] * dIz);
+    TacsScalar dAz = delta_z * dA + w_offset * A;
+    TacsScalar dAy = delta_y * dA;
+    TacsScalar dIy = 0.25 * thickness * width * width +
+                     2.0 * w_offset * delta_z * A + delta_z * delta_z * dA;
+    TacsScalar dIz = 1.0 / 12.0 * t3 + delta_y * delta_y * dA;
+    TacsScalar dIyz = -delta_y * w_offset * A - delta_y * delta_z * dA;
+
+    dfdx[index] += rho * (scale[0] * dA + scale[1] * dAy + scale[2] * dAz +
+                          scale[3] * dIz + scale[4] * dIy - scale[5] * dIyz);
     index++;
   }
   if (thickness_num >= 0) {
     TacsScalar dA = A / thickness;
-    TacsScalar dIy = Iy / thickness;
-    TacsScalar dIz = 3.0 * Iz / thickness;
-
-    dfdx[index] += rho * (scale[0] * dA + scale[3] * dIy + scale[4] * dIz);
+    TacsScalar dAy = delta_y * dA + t_offset * A;
+    TacsScalar dAz = delta_z * dA;
+    TacsScalar dIy = 1.0 / 12.0 * w3 + delta_z * delta_z * dA;
+    TacsScalar dIz = 0.25 * thickness * thickness * width +
+                     2.0 * t_offset * delta_y * A + delta_y * delta_y * dA;
+    TacsScalar dIyz = -delta_z * t_offset * A - delta_y * delta_z * dA;
+
+    dfdx[index] += rho * (scale[0] * dA + scale[1] * dAy + scale[2] * dAz +
+                          scale[3] * dIz + scale[4] * dIy - scale[5] * dIyz);
     index++;
   }
 }
@@ -181,11 +199,14 @@ void TACSIsoRectangleBeamConstitutive::evalStress(int elemIndex,
 
   TacsScalar G = 0.5 * E / (1.0 + nu);
   TacsScalar kcorr = 10.0 * (1.0 + nu) / (12.0 + 11.0 * nu);
+  TacsScalar delta_y = t_offset * thickness;
+  TacsScalar delta_z = w_offset * width;
   TacsScalar A = thickness * width;
   TacsScalar t3 = thickness * thickness * thickness;
   TacsScalar w3 = width * width * width;
-  TacsScalar Iy = 1.0 / 12.0 * thickness * w3;
-  TacsScalar Iz = 1.0 / 12.0 * width * t3;
+  TacsScalar Iy = 1.0 / 12.0 * thickness * w3 + delta_z * delta_z * A;
+  TacsScalar Iz = 1.0 / 12.0 * width * t3 + delta_y * delta_y * A;
+  TacsScalar Iyz = -delta_y * delta_z * A;
   // Torsion constant for rectangle
   TacsScalar a = width / 2.0;
   TacsScalar b = thickness / 2.0;
@@ -194,10 +215,10 @@ void TACSIsoRectangleBeamConstitutive::evalStress(int elemIndex,
   TacsScalar a4 = a * a * a * a;
   TacsScalar J = a * b3 * (16.0 / 3.0 - 3.36 * b / a * (1.0 - b4 / a4 / 12.0));
 
-  s[0] = E * A * e[0];
+  s[0] = E * A * (e[0] - delta_y * e[2] - delta_z * e[3]);
   s[1] = G * J * e[1];
-  s[2] = E * Iz * e[2];
-  s[3] = E * Iy * e[3];
+  s[2] = E * (Iz * e[2] - delta_y * A * e[0] - Iyz * e[3]);
+  s[3] = E * (Iy * e[3] - delta_z * A * e[0] - Iyz * e[2]);
   s[4] = kcorr * G * A * e[4];
   s[5] = kcorr * G * A * e[5];
 }
@@ -209,11 +230,14 @@ void TACSIsoRectangleBeamConstitutive::evalTangentStiffness(
 
   TacsScalar G = 0.5 * E / (1.0 + nu);
   TacsScalar kcorr = 10.0 * (1.0 + nu) / (12.0 + 11.0 * nu);
+  TacsScalar delta_y = t_offset * thickness;
+  TacsScalar delta_z = w_offset * width;
   TacsScalar A = thickness * width;
   TacsScalar t3 = thickness * thickness * thickness;
   TacsScalar w3 = width * width * width;
-  TacsScalar Iy = 1.0 / 12.0 * thickness * w3;
-  TacsScalar Iz = 1.0 / 12.0 * width * t3;
+  TacsScalar Iy = 1.0 / 12.0 * thickness * w3 + delta_z * delta_z * A;
+  TacsScalar Iz = 1.0 / 12.0 * width * t3 + delta_y * delta_y * A;
+  TacsScalar Iyz = -delta_y * delta_z * A;
   // Torsion constant for rectangle
   TacsScalar a = width / 2.0;
   TacsScalar b = thickness / 2.0;
@@ -225,8 +249,11 @@ void TACSIsoRectangleBeamConstitutive::evalTangentStiffness(
   memset(C, 0, NUM_TANGENT_STIFFNESS_ENTRIES * sizeof(TacsScalar));
 
   C[0] = E * A;
+  C[2] = -E * delta_y * A;
+  C[3] = -E * delta_z * A;
   C[6] = G * J;
   C[11] = E * Iz;
+  C[12] = -E * Iyz;
   C[15] = E * Iy;
   C[18] = kcorr * G * A;
   C[20] = kcorr * G * A;
@@ -241,51 +268,70 @@ void TACSIsoRectangleBeamConstitutive::addStressDVSens(
 
   TacsScalar G = 0.5 * E / (1.0 + nu);
   TacsScalar kcorr = 10.0 * (1.0 + nu) / (12.0 + 11.0 * nu);
+  TacsScalar A = thickness * width;
+  TacsScalar delta_y = t_offset * thickness;
+  TacsScalar delta_z = w_offset * width;
+  TacsScalar a = width / 2.0;
+  TacsScalar b = thickness / 2.0;
+  TacsScalar b3 = b * b * b;
+  TacsScalar b4 = b * b3;
+  TacsScalar a4 = a * a * a * a;
 
   int index = 0;
   if (width_num >= 0) {
+    TacsScalar ddelta_z = w_offset;
     TacsScalar dA = thickness;
-    TacsScalar dIz = thickness * thickness * thickness / 12.0;
-    TacsScalar dIy = width * width * thickness / 4.0;
-    TacsScalar a = width / 2.0;
-    TacsScalar da = 1.0 / 2.0;
-    TacsScalar b = thickness / 2.0;
-    TacsScalar b3 = b * b * b;
-    TacsScalar b4 = b * b3;
-    TacsScalar a4 = a * a * a * a;
+    TacsScalar dIz =
+        1.0 / 12.0 * thickness * thickness * thickness + delta_y * delta_y * dA;
+    TacsScalar dIy = width * width * thickness / 4.0 +
+                     2.0 * ddelta_z * delta_z * A + delta_z * delta_z * dA;
+    TacsScalar dIyz = -delta_y * ddelta_z * A + -delta_y * delta_z * dA;
+    TacsScalar da = 0.5;
     TacsScalar da4 = 4.0 * a * a * a * da;
     TacsScalar dJ =
         da * b3 * (16.0 / 3.0 - 3.36 * b / a * (1.0 - b4 / a4 / 12.0)) +
         a * b3 * (3.36 * b / a * da / a * (1.0 - b4 / a4 / 12.0)) +
         a * b3 * (3.36 * b / a * (-b4 / a4 / 12.0 * da4 / a4));
 
-    dfdx[index] += scale * (E * dA * e[0] * psi[0] + G * dJ * e[1] * psi[1] +
-                            E * dIz * e[2] * psi[2] + E * dIy * e[3] * psi[3] +
-                            kcorr * G * dA * e[4] * psi[4] +
-                            kcorr * G * dA * e[5] * psi[5]);
+    dfdx[index] +=
+        scale *
+        (E *
+             (dA * e[0] - delta_y * dA * e[2] -
+              (delta_z * dA + ddelta_z * A) * e[3]) *
+             psi[0] +
+         G * dJ * e[1] * psi[1] +
+         E * (dIz * e[2] - dIyz * e[3] - delta_y * dA * e[0]) * psi[2] +
+         E * (dIy * e[3] - dIyz * e[2] - (delta_z * dA + ddelta_z * A) * e[0]) *
+             psi[3] +
+         kcorr * G * dA * e[4] * psi[4] + kcorr * G * dA * e[5] * psi[5]);
     index++;
   }
   if (thickness_num >= 0) {
+    TacsScalar ddelta_y = t_offset;
     TacsScalar dA = width;
-    TacsScalar dIz = width * thickness * thickness / 4.0;
-    TacsScalar dIy = width * width * width / 12.0;
-    TacsScalar a = width / 2.0;
-    TacsScalar b = thickness / 2.0;
-    TacsScalar db = 1.0 / 2.0;
-    TacsScalar b3 = b * b * b;
+    TacsScalar dIz = width * thickness * thickness / 4.0 +
+                     2.0 * ddelta_y * delta_y * A + delta_y * delta_y * dA;
+    TacsScalar dIy = width * width * width / 12.0 + delta_z * delta_z * dA;
+    TacsScalar dIyz = -ddelta_y * delta_z * A + -delta_y * delta_z * dA;
+    TacsScalar db = 0.5;
     TacsScalar db3 = 3.0 * b * b * db;
-    TacsScalar b4 = b * b3;
     TacsScalar db4 = 4.0 * b3 * db;
-    TacsScalar a4 = a * a * a * a;
     TacsScalar dJ =
         a * db3 * (16.0 / 3.0 - 3.36 * b / a * (1.0 - b4 / a4 / 12.0)) +
         a * b3 * (-3.36 * db / a * (1.0 - b4 / a4 / 12.0)) +
         a * b3 * (-3.36 * b / a * (-db4 / a4 / 12.0));
 
-    dfdx[index] += scale * (E * dA * e[0] * psi[0] + G * dJ * e[1] * psi[1] +
-                            E * dIz * e[2] * psi[2] + E * dIy * e[3] * psi[3] +
-                            kcorr * G * dA * e[4] * psi[4] +
-                            kcorr * G * dA * e[5] * psi[5]);
+    dfdx[index] +=
+        scale *
+        (E *
+             (dA * e[0] - (delta_y * dA + ddelta_y * A) * e[2] -
+              delta_z * dA * e[3]) *
+             psi[0] +
+         G * dJ * e[1] * psi[1] +
+         E * (dIz * e[2] - dIyz * e[3] - (delta_y * dA + ddelta_y * A) * e[0]) *
+             psi[2] +
+         E * (dIy * e[3] - dIyz * e[2] - delta_z * dA * e[0]) * psi[3] +
+         kcorr * G * dA * e[4] * psi[4] + kcorr * G * dA * e[5] * psi[5]);
     index++;
   }
 }
@@ -298,19 +344,20 @@ TacsScalar TACSIsoRectangleBeamConstitutive::evalFailure(int elemIndex,
   TacsScalar e0[6], s0[6];    // 3D stress, NOT beam stresses
   TacsScalar fail_checks[4];  // fail value a four corners
   TacsScalar max_fail = -1e20, ks_sum = 0.0;
-  TacsScalar y_lim[] = {-thickness / 2.0, thickness / 2.0};
-  TacsScalar z_lim[] = {-width / 2.0, width / 2.0};
+  TacsScalar y_lim[] = {-(0.5 + t_offset) * thickness,
+                        (0.5 - t_offset) * thickness};
+  TacsScalar z_lim[] = {-(0.5 + w_offset) * width, (0.5 - w_offset) * width};
 
   int count = 0;
   for (int i = 0; i < 2; i++) {
     for (int j = 0; j < 2; j++, count++) {
       // Compute the combined strain state e0 = [ex, ey, ez, gyz, gxz, gxy]
-      e0[0] = e[0] - y_lim[i] * e[2] - z_lim[j] * e[3];  // ex
+      e0[0] = e[0] + y_lim[i] * e[2] + z_lim[j] * e[3];  // ex
       e0[1] = 0.0;
       e0[2] = 0.0;
       e0[3] = 0.0;
-      e0[4] = e[5] + y_lim[i] * e[1];
-      e0[5] = e[4] - z_lim[j] * e[1];
+      e0[4] = e[5] + (y_lim[i] + t_offset * thickness) * e[1];
+      e0[5] = e[4] - (z_lim[j] + w_offset * width) * e[1];
 
       // Compute the stress
       props->evalStress3D(e0, s0);
@@ -340,8 +387,9 @@ TacsScalar TACSIsoRectangleBeamConstitutive::evalFailureStrainSens(
   TacsScalar fail_checks[4];
   TacsScalar fail_checks_sens[4][6];
   TacsScalar max_fail = -1e20, ks_sum = 0.0;
-  TacsScalar y_lim[] = {-thickness / 2.0, thickness / 2.0};
-  TacsScalar z_lim[] = {-width / 2.0, width / 2.0};
+  TacsScalar y_lim[] = {-(0.5 + t_offset) * thickness,
+                        (0.5 - t_offset) * thickness};
+  TacsScalar z_lim[] = {-(0.5 + w_offset) * width, (0.5 - w_offset) * width};
 
   memset(sens, 0, 6 * sizeof(TacsScalar));
 
@@ -349,12 +397,12 @@ TacsScalar TACSIsoRectangleBeamConstitutive::evalFailureStrainSens(
   for (int i = 0; i < 2; i++) {
     for (int j = 0; j < 2; j++, count++) {
       // Compute the combined strain state e0 = [ex, ey, ez, gyz, gxz, gxy]
-      e0[0] = e[0] - y_lim[i] * e[2] - z_lim[j] * e[3];  // ex
+      e0[0] = e[0] + y_lim[i] * e[2] + z_lim[j] * e[3];  // ex
       e0[1] = 0.0;
       e0[2] = 0.0;
       e0[3] = 0.0;
-      e0[4] = e[5] + y_lim[i] * e[1];
-      e0[5] = e[4] - z_lim[j] * e[1];
+      e0[4] = e[5] + (y_lim[i] + t_offset * thickness) * e[1];
+      e0[5] = e[4] - (z_lim[j] + w_offset * width) * e[1];
 
       // Compute the stress
       props->evalStress3D(e0, s0);
@@ -364,9 +412,11 @@ TacsScalar TACSIsoRectangleBeamConstitutive::evalFailureStrainSens(
 
       props->evalStress3D(dfds0, dfde0);
       fail_checks_sens[count][0] = dfde0[0];
-      fail_checks_sens[count][2] = -y_lim[i] * dfde0[0];
-      fail_checks_sens[count][3] = -z_lim[j] * dfde0[0];
-      fail_checks_sens[count][1] = y_lim[i] * dfde0[4] - z_lim[j] * dfde0[5];
+      fail_checks_sens[count][2] = y_lim[i] * dfde0[0];
+      fail_checks_sens[count][3] = z_lim[j] * dfde0[0];
+      fail_checks_sens[count][1] =
+          (y_lim[i] + t_offset * thickness) * dfde0[4] -
+          (z_lim[j] + w_offset * width) * dfde0[5];
       fail_checks_sens[count][5] = dfde0[4];
       fail_checks_sens[count][4] = dfde0[5];
 
@@ -405,30 +455,33 @@ void TACSIsoRectangleBeamConstitutive::addFailureDVSens(
     TacsScalar fail_checks[4];
     TacsScalar fail_checks_sens[4];
     TacsScalar max_fail = -1e20, ks_sum = 0.0;
-    TacsScalar y_lim[] = {-thickness / 2.0, thickness / 2.0};
-    TacsScalar z_lim[] = {-width / 2.0, width / 2.0};
-    TacsScalar dy_lim[] = {0.0, 0.0};
-    TacsScalar dz_lim[] = {0.0, 0.0};
+    TacsScalar y_lim[] = {-(0.5 + t_offset) * thickness,
+                          (0.5 - t_offset) * thickness};
+    TacsScalar z_lim[] = {-(0.5 + w_offset) * width, (0.5 - w_offset) * width};
+    TacsScalar dwidth = 0.0, dthickness = 0.0;
 
     if (dvNum == width_num) {
-      dz_lim[0] = -0.5;
-      dz_lim[1] = 0.5;
+      dwidth = 1.0;
     }
     if (dvNum == thickness_num) {
-      dy_lim[0] = -0.5;
-      dy_lim[1] = 0.5;
+      dthickness = 1.0;
     }
 
+    TacsScalar dy_lim[] = {-(0.5 + t_offset) * dthickness,
+                           (0.5 - t_offset) * dthickness};
+    TacsScalar dz_lim[] = {-(0.5 + w_offset) * dwidth,
+                           (0.5 - w_offset) * dwidth};
+
     int count = 0;
     for (int i = 0; i < 2; i++) {
       for (int j = 0; j < 2; j++, count++) {
         // Compute the combined strain state e0 = [ex, ey, ez, gyz, gxz, gxy]
-        e0[0] = e[0] - y_lim[i] * e[2] - z_lim[j] * e[3];  // ex
+        e0[0] = e[0] + y_lim[i] * e[2] + z_lim[j] * e[3];  // ex
         e0[1] = 0.0;
         e0[2] = 0.0;
         e0[3] = 0.0;
-        e0[4] = e[5] + y_lim[i] * e[1];
-        e0[5] = e[4] - z_lim[j] * e[1];
+        e0[4] = e[5] + (y_lim[i] + t_offset * thickness) * e[1];
+        e0[5] = e[4] - (z_lim[j] + w_offset * width) * e[1];
 
         // Compute the stress
         props->evalStress3D(e0, s0);
@@ -438,8 +491,10 @@ void TACSIsoRectangleBeamConstitutive::addFailureDVSens(
         props->evalStress3D(s0d, e0d);
 
         fail_checks_sens[count] =
-            e0d[0] * (-dy_lim[i] * e[2] - dz_lim[j] * e[3]) +
-            (e0d[4] * dy_lim[i] - e0d[5] * dz_lim[j]) * e[1];
+            e0d[0] * (dy_lim[i] * e[2] + dz_lim[j] * e[3]) +
+            (e0d[4] * (dy_lim[i] + t_offset * dthickness) -
+             e0d[5] * (dz_lim[j] + w_offset * dwidth)) *
+                e[1];
 
         if (TacsRealPart(fail_checks[count]) > TacsRealPart(max_fail)) {
           max_fail = fail_checks[count];
diff --git a/src/constitutive/TACSIsoRectangleBeamConstitutive.h b/src/constitutive/TACSIsoRectangleBeamConstitutive.h
index 8fdc8b978..24f814500 100644
--- a/src/constitutive/TACSIsoRectangleBeamConstitutive.h
+++ b/src/constitutive/TACSIsoRectangleBeamConstitutive.h
@@ -40,7 +40,9 @@ class TACSIsoRectangleBeamConstitutive : public TACSBeamConstitutive {
                                    int _width_num, int _thickness_num,
                                    TacsScalar _lb_width, TacsScalar _ub_width,
                                    TacsScalar _lb_thickness,
-                                   TacsScalar _ub_thickness);
+                                   TacsScalar _ub_thickness,
+                                   TacsScalar _w_offset = 0.0,
+                                   TacsScalar _t_offset = 0.0);
   ~TACSIsoRectangleBeamConstitutive();
 
   // Retrieve the global design variable numbers
@@ -118,6 +120,7 @@ class TACSIsoRectangleBeamConstitutive : public TACSBeamConstitutive {
   TacsScalar lb_thickness, ub_thickness;
   TacsScalar lb_width, ub_width;
   TacsScalar ks_weight;
+  TacsScalar w_offset, t_offset;
   // The object name
   static const char *constName;
 };
diff --git a/src/constitutive/TACSIsoShellConstitutive.cpp b/src/constitutive/TACSIsoShellConstitutive.cpp
index 499b5ed2e..3fc0a8395 100644
--- a/src/constitutive/TACSIsoShellConstitutive.cpp
+++ b/src/constitutive/TACSIsoShellConstitutive.cpp
@@ -27,7 +27,7 @@ const char *TACSIsoShellConstitutive::constName = "TACSIsoShellConstitutive";
 */
 TACSIsoShellConstitutive::TACSIsoShellConstitutive(
     TACSMaterialProperties *props, TacsScalar _t, int _tNum, TacsScalar _tlb,
-    TacsScalar _tub) {
+    TacsScalar _tub, TacsScalar _tOffset) {
   properties = props;
   if (properties) {
     properties->incref();
@@ -37,6 +37,7 @@ TACSIsoShellConstitutive::TACSIsoShellConstitutive(
   tNum = _tNum;
   tlb = _tlb;
   tub = _tub;
+  tOffset = _tOffset;
   kcorr = 5.0 / 6.0;
   ksWeight = 100.0;
 }
@@ -122,8 +123,8 @@ void TACSIsoShellConstitutive::evalMassMoments(int elemIndex, const double pt[],
   if (properties) {
     TacsScalar rho = properties->getDensity();
     moments[0] = rho * t;
-    moments[1] = 0.0;
-    moments[2] = rho * t * t * t / 12.0;
+    moments[1] = -rho * t * t * tOffset;
+    moments[2] = rho * t * t * t * (tOffset * tOffset + 1.0 / 12.0);
   }
 }
 
@@ -133,7 +134,8 @@ void TACSIsoShellConstitutive::addMassMomentsDVSens(
     const TacsScalar scale[], int dvLen, TacsScalar dfdx[]) {
   if (properties && tNum >= 0) {
     TacsScalar rho = properties->getDensity();
-    dfdx[0] += rho * (scale[0] + 0.25 * t * t * scale[2]);
+    dfdx[0] += rho * (scale[0] - 2.0 * tOffset * t * scale[1] +
+                      (3.0 * tOffset * tOffset + 0.25) * t * t * scale[2]);
   }
 }
 
@@ -147,7 +149,7 @@ TacsScalar TACSIsoShellConstitutive::evalSpecificHeat(int elemIndex,
   return 0.0;
 }
 
-// Evaluate the stresss
+// Evaluate the stress
 void TACSIsoShellConstitutive::evalStress(int elemIndex, const double pt[],
                                           const TacsScalar X[],
                                           const TacsScalar e[],
@@ -167,6 +169,8 @@ void TACSIsoShellConstitutive::evalStress(int elemIndex, const double pt[],
     for (int i = 0; i < 6; i++) {
       D[i] = I * A[i];
       A[i] *= t;
+      B[i] += -tOffset * t * A[i];
+      D[i] += tOffset * tOffset * t * t * A[i];
     }
 
     // Set the through-thickness shear stiffness
@@ -207,6 +211,8 @@ void TACSIsoShellConstitutive::evalTangentStiffness(int elemIndex,
     for (int i = 0; i < 6; i++) {
       D[i] = I * A[i];
       A[i] *= t;
+      B[i] += -tOffset * t * A[i];
+      D[i] += tOffset * tOffset * t * t * A[i];
     }
 
     // Set the through-thickness shear stiffness
@@ -231,13 +237,16 @@ void TACSIsoShellConstitutive::addStressDVSens(int elemIndex, TacsScalar scale,
     TacsScalar A[6];
     properties->evalTangentStiffness2D(A);
 
-    TacsScalar dI = t * t / 4.0;
+    TacsScalar dI = (3.0 * tOffset * tOffset + 0.25) * t * t;
 
-    dfdx[0] += scale * (mat3x3SymmInner(A, &psi[0], &e[0]) +
-                        dI * mat3x3SymmInner(A, &psi[3], &e[3]) +
-                        (5.0 / 6.0) * A[5] *
-                            (psi[6] * e[6] + psi[7] * e[7] +
-                             DRILLING_REGULARIZATION * psi[8] * e[8]));
+    dfdx[0] +=
+        scale * (mat3x3SymmInner(A, &psi[0], &e[0]) +
+                 dI * mat3x3SymmInner(A, &psi[3], &e[3]) +
+                 -2.0 * tOffset * t * mat3x3SymmInner(A, &psi[0], &e[3]) +
+                 -2.0 * tOffset * t * mat3x3SymmInner(A, &psi[3], &e[0]) +
+                 (5.0 / 6.0) * A[5] *
+                     (psi[6] * e[6] + psi[7] * e[7] +
+                      DRILLING_REGULARIZATION * psi[8] * e[8]));
   }
 }
 
@@ -248,15 +257,16 @@ TacsScalar TACSIsoShellConstitutive::evalFailure(int elemIndex,
                                                  const TacsScalar e[]) {
   if (properties) {
     TacsScalar et[3], eb[3];
-    TacsScalar ht = 0.5 * t;
+    TacsScalar ht = (0.5 - tOffset) * t;
+    TacsScalar hb = -(0.5 + tOffset) * t;
 
     et[0] = e[0] + ht * e[3];
     et[1] = e[1] + ht * e[4];
     et[2] = e[2] + ht * e[5];
 
-    eb[0] = e[0] - ht * e[3];
-    eb[1] = e[1] - ht * e[4];
-    eb[2] = e[2] - ht * e[5];
+    eb[0] = e[0] + hb * e[3];
+    eb[1] = e[1] + hb * e[4];
+    eb[2] = e[2] + hb * e[5];
 
     TacsScalar C[6];
     properties->evalTangentStiffness2D(C);
@@ -298,15 +308,16 @@ TacsScalar TACSIsoShellConstitutive::evalFailureStrainSens(int elemIndex,
 
   if (properties) {
     TacsScalar et[3], eb[3];
-    TacsScalar ht = 0.5 * t;
+    TacsScalar ht = (0.5 - tOffset) * t;
+    TacsScalar hb = -(0.5 + tOffset) * t;
 
     et[0] = e[0] + ht * e[3];
     et[1] = e[1] + ht * e[4];
     et[2] = e[2] + ht * e[5];
 
-    eb[0] = e[0] - ht * e[3];
-    eb[1] = e[1] - ht * e[4];
-    eb[2] = e[2] - ht * e[5];
+    eb[0] = e[0] + hb * e[3];
+    eb[1] = e[1] + hb * e[4];
+    eb[2] = e[2] + hb * e[5];
 
     TacsScalar C[6];
     properties->evalTangentStiffness2D(C);
@@ -354,9 +365,9 @@ TacsScalar TACSIsoShellConstitutive::evalFailureStrainSens(int elemIndex,
     sens[1] += ksFactor * phi[1];
     sens[2] += ksFactor * phi[2];
 
-    sens[3] -= ksFactor * ht * phi[0];
-    sens[4] -= ksFactor * ht * phi[1];
-    sens[5] -= ksFactor * ht * phi[2];
+    sens[3] += ksFactor * hb * phi[0];
+    sens[4] += ksFactor * hb * phi[1];
+    sens[5] += ksFactor * hb * phi[2];
 
     sens[6] = sens[7] = sens[8] = 0.0;
 
@@ -374,15 +385,16 @@ void TACSIsoShellConstitutive::addFailureDVSens(int elemIndex, TacsScalar scale,
                                                 TacsScalar dfdx[]) {
   if (properties && tNum >= 0 && dvLen >= 1) {
     TacsScalar et[3], eb[3];
-    TacsScalar ht = 0.5 * t;
+    TacsScalar ht = (0.5 - tOffset) * t;
+    TacsScalar hb = -(0.5 + tOffset) * t;
 
     et[0] = e[0] + ht * e[3];
     et[1] = e[1] + ht * e[4];
     et[2] = e[2] + ht * e[5];
 
-    eb[0] = e[0] - ht * e[3];
-    eb[1] = e[1] - ht * e[4];
-    eb[2] = e[2] - ht * e[5];
+    eb[0] = e[0] + hb * e[3];
+    eb[1] = e[1] + hb * e[4];
+    eb[2] = e[2] + hb * e[5];
 
     TacsScalar C[6];
     properties->evalTangentStiffness2D(C);
@@ -410,14 +422,14 @@ void TACSIsoShellConstitutive::addFailureDVSens(int elemIndex, TacsScalar scale,
     properties->vonMisesFailure2DStressSens(st, psi);
     mat3x3SymmMult(C, psi, phi);
     TacsScalar ksFactor = exp(ksWeight * (top - ksMax)) / ksSum;
-    dfdx[0] += ksFactor * 0.5 * scale *
+    dfdx[0] += ksFactor * (0.5 - tOffset) * scale *
                (phi[0] * e[3] + phi[1] * e[4] + phi[2] * e[5]);
 
     // Contribution from plate bottom
     properties->vonMisesFailure2DStressSens(sb, psi);
     mat3x3SymmMult(C, psi, phi);
     ksFactor = exp(ksWeight * (bottom - ksMax)) / ksSum;
-    dfdx[0] -= ksFactor * 0.5 * scale *
+    dfdx[0] += ksFactor * -(0.5 + tOffset) * scale *
                (phi[0] * e[3] + phi[1] * e[4] + phi[2] * e[5]);
   }
 }
diff --git a/src/constitutive/TACSIsoShellConstitutive.h b/src/constitutive/TACSIsoShellConstitutive.h
index 27ac4e4fa..f38386471 100644
--- a/src/constitutive/TACSIsoShellConstitutive.h
+++ b/src/constitutive/TACSIsoShellConstitutive.h
@@ -32,7 +32,7 @@ class TACSIsoShellConstitutive : public TACSShellConstitutive {
  public:
   TACSIsoShellConstitutive(TACSMaterialProperties *props, TacsScalar _t = 1.0,
                            int _tNum = -1, TacsScalar _tlb = 0.0,
-                           TacsScalar _tub = 1.0);
+                           TacsScalar _tub = 1.0, TacsScalar _tOffset = 0.0);
   ~TACSIsoShellConstitutive();
 
   // Retrieve the global design variable numbers
@@ -126,7 +126,7 @@ class TACSIsoShellConstitutive : public TACSShellConstitutive {
 
   // Store information about the design variable
   TacsScalar kcorr;  // The shear correction factor
-  TacsScalar t, tlb, tub;
+  TacsScalar t, tlb, tub, tOffset;
   int tNum;
   TacsScalar ksWeight;  // ks weight used in failure calc
 
diff --git a/src/constitutive/TACSSmearedCompositeShellConstitutive.cpp b/src/constitutive/TACSSmearedCompositeShellConstitutive.cpp
index 80f983f22..30f7c8f50 100644
--- a/src/constitutive/TACSSmearedCompositeShellConstitutive.cpp
+++ b/src/constitutive/TACSSmearedCompositeShellConstitutive.cpp
@@ -32,7 +32,8 @@ TACSSmearedCompositeShellConstitutive::TACSSmearedCompositeShellConstitutive(
     const TacsScalar *_ply_angles, const TacsScalar *_ply_fractions,
     int _thickness_dv_num, const int *_ply_fraction_dv_nums,
     TacsScalar _thickness_lb, TacsScalar _thickness_ub,
-    const TacsScalar *_ply_fraction_lb, const TacsScalar *_ply_fraction_ub) {
+    const TacsScalar *_ply_fraction_lb, const TacsScalar *_ply_fraction_ub,
+    TacsScalar _t_offset) {
   num_plies = _num_plies;
 
   thickness = _thickness;
@@ -76,6 +77,7 @@ TACSSmearedCompositeShellConstitutive::TACSSmearedCompositeShellConstitutive(
   }
 
   ks_weight = 100.0;
+  t_offset = _t_offset;
   nfvals = 2 * num_plies;
   fvals = new TacsScalar[nfvals];
   dks_vals = new TacsScalar[nfvals];
@@ -236,12 +238,15 @@ void TACSSmearedCompositeShellConstitutive::evalMassMoments(
   moments[2] = 0.0;
 
   // Compute the contribution to the mass moment from each layer
+  TacsScalar t2 = thickness * thickness;
+  TacsScalar t3 = t2 * thickness;
   for (int i = 0; i < num_plies; i++) {
     TacsScalar rho_ply = ply_props[i]->getDensity();
 
     moments[0] += thickness * rho_ply * ply_fractions[i];
+    moments[1] -= t_offset * t2 * rho_ply * ply_fractions[i];
     moments[2] +=
-        thickness * thickness * thickness * rho_ply * ply_fractions[i] / 12.0;
+        (t_offset * t_offset + 1.0 / 12.0) * t3 * rho_ply * ply_fractions[i];
   }
 }
 
@@ -254,8 +259,10 @@ void TACSSmearedCompositeShellConstitutive::addMassMomentsDVSens(
     for (int i = 0; i < num_plies; i++) {
       TacsScalar rho_ply = ply_props[i]->getDensity();
 
-      dfdx[index] +=
-          rho_ply * (scale[0] + 0.25 * thickness * thickness * scale[2]);
+      dfdx[index] += rho_ply * ply_fractions[i] *
+                     (scale[0] - 2.0 * t_offset * thickness * scale[1] +
+                      (3.0 * t_offset * t_offset + 0.25) * thickness *
+                          thickness * scale[2]);
     }
     index++;
   }
@@ -263,9 +270,10 @@ void TACSSmearedCompositeShellConstitutive::addMassMomentsDVSens(
     if (ply_fraction_dv_nums[i] >= 0) {
       TacsScalar rho_ply = ply_props[i]->getDensity();
 
-      dfdx[index] +=
-          rho_ply * (scale[0] * thickness +
-                     thickness * thickness * thickness * scale[2] / 12.0);
+      dfdx[index] += rho_ply * (scale[0] * thickness -
+                                t_offset * thickness * thickness * scale[1] +
+                                (t_offset * t_offset + 1.0 / 12.0) * thickness *
+                                    thickness * thickness * scale[2]);
       index++;
     }
   }
@@ -299,11 +307,13 @@ TacsScalar TACSSmearedCompositeShellConstitutive::evalFSDTStiffness(
     ply_props[k]->calculateAbar(ply_angles[k], Abar);
 
     TacsScalar a = ply_fractions[k] * thickness;
-    TacsScalar d =
-        1.0 / 12.0 * ply_fractions[k] * (thickness * thickness * thickness);
+    TacsScalar b = t_offset * ply_fractions[k] * thickness * thickness;
+    TacsScalar d = (t_offset * t_offset + 1.0 / 12.0) * ply_fractions[k] *
+                   (thickness * thickness * thickness);
 
     for (int i = 0; i < 6; i++) {
       A[i] += a * Qbar[i];
+      B[i] += -b * Qbar[i];
       D[i] += d * Qbar[i];
     }
 
@@ -335,7 +345,8 @@ void TACSSmearedCompositeShellConstitutive::addStressDVSens(
     int elemIndex, TacsScalar scale, const double pt[], const TacsScalar X[],
     const TacsScalar e[], const TacsScalar psi[], int dvLen,
     TacsScalar dfdx[]) {
-  TacsScalar dA[6] = {0.0}, dD[6] = {0.0}, dAs[3] = {0.0}, ddrill;
+  TacsScalar dA[6] = {0.0}, dB[6] = {0.0}, dD[6] = {0.0}, dAs[3] = {0.0},
+             ddrill;
 
   int index = 0;
   if (thickness_dv_num >= 0) {
@@ -345,10 +356,13 @@ void TACSSmearedCompositeShellConstitutive::addStressDVSens(
       ply_props[k]->calculateAbar(ply_angles[k], Abar);
 
       TacsScalar da = ply_fractions[k];
-      TacsScalar dd = 0.25 * ply_fractions[k] * (thickness * thickness);
+      TacsScalar db = 2.0 * t_offset * ply_fractions[k] * thickness;
+      TacsScalar dd = (3.0 * t_offset * t_offset + 0.25) * ply_fractions[k] *
+                      (thickness * thickness);
 
       for (int i = 0; i < 6; i++) {
         dA[i] += da * Qbar[i];
+        dB[i] += -db * Qbar[i];
         dD[i] += dd * Qbar[i];
       }
 
@@ -362,6 +376,8 @@ void TACSSmearedCompositeShellConstitutive::addStressDVSens(
     dfdx[index] +=
         scale * (mat3x3SymmInner(dA, &psi[0], &e[0]) +
                  mat3x3SymmInner(dD, &psi[3], &e[3]) +
+                 mat3x3SymmInner(dB, &psi[0], &e[3]) +
+                 mat3x3SymmInner(dB, &psi[3], &e[0]) +
                  mat2x2SymmInner(dAs, &psi[6], &e[6]) + ddrill * psi[8] * e[8]);
     index++;
   }
@@ -373,10 +389,13 @@ void TACSSmearedCompositeShellConstitutive::addStressDVSens(
       ply_props[k]->calculateAbar(ply_angles[k], Abar);
 
       TacsScalar da = thickness;
-      TacsScalar dd = 1.0 / 12.0 * (thickness * thickness * thickness);
+      TacsScalar db = t_offset * thickness * thickness;
+      TacsScalar dd = (t_offset * t_offset + 1.0 / 12.0) *
+                      (thickness * thickness * thickness);
 
       for (int i = 0; i < 6; i++) {
         dA[i] = da * Qbar[i];
+        dB[i] = -db * Qbar[i];
         dD[i] = dd * Qbar[i];
       }
 
@@ -388,6 +407,8 @@ void TACSSmearedCompositeShellConstitutive::addStressDVSens(
 
       dfdx[index] += scale * (mat3x3SymmInner(dA, &psi[0], &e[0]) +
                               mat3x3SymmInner(dD, &psi[3], &e[3]) +
+                              mat3x3SymmInner(dB, &psi[0], &e[3]) +
+                              mat3x3SymmInner(dB, &psi[3], &e[0]) +
                               mat2x2SymmInner(dAs, &psi[6], &e[6]) +
                               ddrill * psi[8] * e[8]);
       index++;
@@ -413,8 +434,8 @@ TacsScalar TACSSmearedCompositeShellConstitutive::evalFailure(
 void TACSSmearedCompositeShellConstitutive::evalPlyTopBottomFailure(
     const TacsScalar strain[], TacsScalar fvals[]) {
   // Compute the total thickness of the laminate
-  TacsScalar tb = -0.5 * thickness;
-  TacsScalar tt = 0.5 * thickness;
+  TacsScalar tb = (-0.5 - t_offset) * thickness;
+  TacsScalar tt = (0.5 - t_offset) * thickness;
 
   for (int i = 0; i < num_plies; i++) {
     TacsScalar lamStrain[3] = {0.0};
@@ -436,8 +457,8 @@ TacsScalar TACSSmearedCompositeShellConstitutive::evalFailureStrainSens(
   sens[6] = sens[7] = sens[8] = 0.0;
 
   // Compute the total thickness of the laminate
-  TacsScalar tb = -0.5 * thickness;
-  TacsScalar tt = 0.5 * thickness;
+  TacsScalar tb = (-0.5 - t_offset) * thickness;
+  TacsScalar tt = (0.5 - t_offset) * thickness;
 
   evalPlyTopBottomFailure(strain, fvals);
 
@@ -482,8 +503,8 @@ void TACSSmearedCompositeShellConstitutive::addFailureDVSens(
   int index = 0;
   if (thickness_dv_num >= 0 && dvLen >= 1) {
     // Compute the total thickness of the laminate
-    TacsScalar tb = -0.5 * thickness;
-    TacsScalar tt = 0.5 * thickness;
+    TacsScalar tb = (-0.5 - t_offset) * thickness;
+    TacsScalar tt = (0.5 - t_offset) * thickness;
 
     evalPlyTopBottomFailure(e, fvals);
     TacsScalar ks_val = ksAggregationSens(fvals, nfvals, ks_weight, dks_vals);
@@ -495,14 +516,14 @@ void TACSSmearedCompositeShellConstitutive::addFailureDVSens(
       getLaminaStrain(e, tb, lamStrain);
       TacsScalar fval_bot =
           ply_props[i]->failureStrainSens(ply_angles[i], lamStrain, phi);
-      dfdx[index] += -0.5 * scale * dks_vals[2 * i + 0] *
+      dfdx[index] += (-0.5 - t_offset) * scale * dks_vals[2 * i + 0] *
                      (phi[0] * e[3] + phi[1] * e[4] + phi[2] * e[5]);
 
       // plate top stress
       getLaminaStrain(e, tt, lamStrain);
       TacsScalar fval_top =
           ply_props[i]->failureStrainSens(ply_angles[i], lamStrain, phi);
-      dfdx[index] += 0.5 * scale * dks_vals[2 * i + 1] *
+      dfdx[index] += (0.5 - t_offset) * scale * dks_vals[2 * i + 1] *
                      (phi[0] * e[3] + phi[1] * e[4] + phi[2] * e[5]);
     }
 
@@ -587,3 +608,10 @@ void TACSSmearedCompositeShellConstitutive::getPlyFractions(
     _ply_fractions[i] = ply_fractions[i];
   }
 }
+
+/*
+  Get thickness offset
+*/
+TacsScalar TACSSmearedCompositeShellConstitutive::getThicknessOffset() {
+  return t_offset;
+}
diff --git a/src/constitutive/TACSSmearedCompositeShellConstitutive.h b/src/constitutive/TACSSmearedCompositeShellConstitutive.h
index e5a62ddd9..f9cb084d1 100644
--- a/src/constitutive/TACSSmearedCompositeShellConstitutive.h
+++ b/src/constitutive/TACSSmearedCompositeShellConstitutive.h
@@ -30,7 +30,7 @@ class TACSSmearedCompositeShellConstitutive : public TACSShellConstitutive {
       int _thickness_dv_num = -1, const int *_ply_fraction_dv_nums = NULL,
       TacsScalar _thickness_lb = 0.0, TacsScalar _thickness_ub = 1e20,
       const TacsScalar *_ply_fraction_lb = NULL,
-      const TacsScalar *_ply_fraction_ub = NULL);
+      const TacsScalar *_ply_fraction_ub = NULL, TacsScalar _t_offset = 0.0);
   ~TACSSmearedCompositeShellConstitutive();
 
   // Retrieve the global design variable numbers
@@ -117,6 +117,7 @@ class TACSSmearedCompositeShellConstitutive : public TACSShellConstitutive {
   TacsScalar getLaminateThickness();
   void getPlyAngles(TacsScalar *_ply_angles);
   void getPlyFractions(TacsScalar *_ply_fractions);
+  TacsScalar getThicknessOffset();
 
  private:
   // Store information about the design variable
@@ -133,6 +134,7 @@ class TACSSmearedCompositeShellConstitutive : public TACSShellConstitutive {
   int nfvals;
   TacsScalar *fvals, *dks_vals;
   TacsScalar **dfvals;
+  TacsScalar t_offset;
 
   // The object name
   static const char *constName;
diff --git a/src/elements/shell/TACSBeamCentrifugalForce.h b/src/elements/shell/TACSBeamCentrifugalForce.h
index 96a36e6f8..7e095208e 100644
--- a/src/elements/shell/TACSBeamCentrifugalForce.h
+++ b/src/elements/shell/TACSBeamCentrifugalForce.h
@@ -4,6 +4,7 @@
 #include "TACSBeamConstitutive.h"
 #include "TACSBeamElementBasis.h"
 #include "TACSBeamElementQuadrature.h"
+#include "TACSBeamElementTransform.h"
 #include "TACSBeamUtilities.h"
 #include "TACSElement.h"
 #include "TACSElementAlgebra.h"
@@ -14,9 +15,12 @@
 template <int vars_per_node, class quadrature, class basis>
 class TACSBeamCentrifugalForce : public TACSElement {
  public:
-  TACSBeamCentrifugalForce(TACSBeamConstitutive *_con,
+  TACSBeamCentrifugalForce(TACSBeamTransform *_transform,
+                           TACSBeamConstitutive *_con,
                            const TacsScalar _omegaVec[],
                            const TacsScalar _rotCenter[]) {
+    transform = _transform;
+    transform->incref();
     con = _con;
     con->incref();
     memcpy(omegaVec, _omegaVec, 3 * sizeof(TacsScalar));
@@ -24,6 +28,9 @@ class TACSBeamCentrifugalForce : public TACSElement {
   }
 
   ~TACSBeamCentrifugalForce() {
+    if (transform) {
+      transform->decref();
+    }
     if (con) {
       con->decref();
     }
@@ -80,6 +87,13 @@ class TACSBeamCentrifugalForce : public TACSElement {
     // Compute the number of quadrature points
     const int nquad = quadrature::getNumQuadraturePoints();
 
+    // Get the reference axis
+    const A2D::Vec3 &axis = transform->getRefAxis();
+
+    // Compute the normal directions
+    TacsScalar fn1[3 * basis::NUM_NODES], fn2[3 * basis::NUM_NODES];
+    TacsBeamComputeNodeNormals<basis>(Xpts, axis, fn1, fn2);
+
     // Loop over each quadrature point and add the residual contribution
     for (int quad_index = 0; quad_index < nquad; quad_index++) {
       // Get the quadrature weight
@@ -89,15 +103,22 @@ class TACSBeamCentrifugalForce : public TACSElement {
       // Tangent to the beam
       A2D::Vec3 X0, X0xi;
 
+      // Interpolated normal directions
+      A2D::ADVec3 n1, n2;
+
       // Compute X, X,xi and the interpolated normal
       basis::template interpFields<3, 3>(pt, Xpts, X0.x);
       basis::template interpFieldsGrad<3, 3>(pt, Xpts, X0xi.x);
+      basis::template interpFields<3, 3>(pt, fn1, n1.x);
+      basis::template interpFields<3, 3>(pt, fn2, n2.x);
 
       // Compute the determinant of the transform
       A2D::Scalar detXd;
       A2D::Vec3Norm(X0xi, detXd);
 
-      TacsScalar mass = con->evalDensity(elemIndex, pt, X0.x);
+      TacsScalar mass_moment[6];
+      con->evalMassMoments(elemIndex, pt, X0.x, mass_moment);
+      TacsScalar mass = mass_moment[0];
 
       TacsScalar r[3], wxr[3], ac[3];
 
@@ -106,6 +127,11 @@ class TACSBeamCentrifugalForce : public TACSElement {
       r[1] = X0.x[1] - rotCenter[1];
       r[2] = X0.x[2] - rotCenter[2];
 
+      // Account for beam offset
+      for (int i = 0; i < 3; i++) {
+        r[i] -= (mass_moment[1] * n1.x[i] + mass_moment[2] * n2.x[i]) / mass;
+      }
+
       // Compute omega x r
       crossProduct(omegaVec, r, wxr);
 
@@ -113,17 +139,22 @@ class TACSBeamCentrifugalForce : public TACSElement {
       crossProduct(omegaVec, wxr, ac);
 
       // Compute the traction
-      TacsScalar tr[3];
+      TacsScalar tr[vars_per_node] = {0.0};
       tr[0] = detXd.value * weight * mass * ac[0];
       tr[1] = detXd.value * weight * mass * ac[1];
       tr[2] = detXd.value * weight * mass * ac[2];
+      // Add moment terms if theres a beam offset
+      crossProductAdd(-detXd.value * weight * mass_moment[1], n1.x, ac, &tr[3]);
+      crossProductAdd(-detXd.value * weight * mass_moment[2], n2.x, ac, &tr[3]);
 
-      basis::template addInterpFieldsTranspose<vars_per_node, 3>(pt, tr, res);
+      basis::template addInterpFieldsTranspose<vars_per_node, vars_per_node>(
+          pt, tr, res);
     }
   }
 
  private:
   TacsScalar omegaVec[3], rotCenter[3];
+  TACSBeamTransform *transform;
   TACSBeamConstitutive *con;
 };
 
diff --git a/src/elements/shell/TACSBeamElement.h b/src/elements/shell/TACSBeamElement.h
index 95759e2a4..26cd69a6e 100644
--- a/src/elements/shell/TACSBeamElement.h
+++ b/src/elements/shell/TACSBeamElement.h
@@ -6,6 +6,7 @@
 #include "TACSBeamElementBasis.h"
 #include "TACSBeamElementModel.h"
 #include "TACSBeamElementQuadrature.h"
+#include "TACSBeamElementTransform.h"
 #include "TACSBeamInertialForce.h"
 #include "TACSBeamTraction.h"
 #include "TACSBeamUtilities.h"
@@ -15,120 +16,6 @@
 #include "TACSGaussQuadrature.h"
 #include "a2d.h"
 
-/*
-  Compute the transformation from the local coordinates
-*/
-class TACSBeamTransform : public TACSObject {
- public:
-  /*
-    Given the local beam element reference frame Xf, compute the
-    transformation from the global coordinates to the shell-aligned local axis.
-  */
-  virtual void computeTransform(const TacsScalar Xxi[], TacsScalar T[]) = 0;
-  virtual void addTransformSens(const TacsScalar X0xi_vals[],
-                                const TacsScalar dTvals[],
-                                TacsScalar dX0xi[]) = 0;
-  virtual A2D::Vec3 &getRefAxis() = 0;
-};
-
-/*
-  Compute the transformation
-*/
-class TACSBeamRefAxisTransform : public TACSBeamTransform {
- public:
-  TACSBeamRefAxisTransform(const TacsScalar axis_dir[]) {
-    A2D::Vec3 axdir(axis_dir);
-    A2D::Vec3Normalize normalize(axdir, axis);
-  }
-
-  void computeTransform(const TacsScalar X0xi_vals[], TacsScalar Tvals[]) {
-    // Normalize the first direction.
-    A2D::Vec3 X0xi(X0xi_vals);
-    A2D::Vec3 t1;
-    A2D::Vec3Normalize normalizet1(X0xi, t1);
-
-    // t2_dir = axis - dot(t1, axis) * t1
-    A2D::Vec3 t2_dir;
-    A2D::Scalar dot;
-    A2D::Vec3Dot dott1(axis, t1, dot);
-    A2D::Vec3Axpy axpy(-1.0, dot, t1, axis, t2_dir);
-
-    // Check if ref axis is parallel to beam
-    if (abs(TacsRealPart(dot.value)) > 1.0 - SMALL_NUM) {
-      fprintf(stderr,
-              "TACSBeamRefAxisTransform: Error, user-provided reference axis "
-              "is parallel to beam axis. "
-              "Element behavior may be ill-conditioned.\n");
-    }
-
-    // Compute the t2 direction
-    A2D::Vec3 t2;
-    A2D::Vec3Normalize normalizet2(t2_dir, t2);
-
-    // Compute the n2 direction
-    A2D::Vec3 t3;
-    A2D::Vec3CrossProduct cross(t1, t2, t3);
-
-    // Assemble the reference frame
-    A2D::Mat3x3 T;
-    A2D::Mat3x3FromThreeVec3 assembleT(t1, t2, t3, T);
-
-    for (int i = 0; i < 9; i++) {
-      Tvals[i] = T.A[i];
-    }
-  }
-
-  void addTransformSens(const TacsScalar X0xi_vals[], const TacsScalar dTvals[],
-                        TacsScalar dX0xi[]) {
-    // Normalize the first direction.
-    A2D::ADVec3 X0xi(X0xi_vals);
-    A2D::ADVec3 t1;
-    A2D::ADVec3Normalize normalizet1(X0xi, t1);
-
-    // t2_dir = axis - dot(t1, axis) * t1
-    A2D::ADVec3 t2_dir;
-    A2D::ADScalar dot;
-    A2D::Vec3ADVecDot dott1(axis, t1, dot);
-    A2D::ADVec3VecADScalarAxpy axpy(-1.0, dot, t1, axis, t2_dir);
-
-    // Compute the t2 direction
-    A2D::ADVec3 t2;
-    A2D::ADVec3Normalize normalizet2(t2_dir, t2);
-
-    // Compute the n2 direction
-    A2D::ADVec3 t3;
-    A2D::ADVec3CrossProduct cross(t1, t2, t3);
-
-    // Assemble the referece frame
-    A2D::ADMat3x3 T(NULL, dTvals);  // Set the seeds for T
-    A2D::ADMat3x3FromThreeADVec3 assembleT(t1, t2, t3, T);
-
-    // Reverse the operations to get the derivative w.r.t. X0
-    assembleT.reverse();
-    cross.reverse();
-    normalizet2.reverse();
-    axpy.reverse();
-    dott1.reverse();
-    normalizet1.reverse();
-
-    for (int i = 0; i < 3; i++) {
-      dX0xi[i] += X0xi.xd[i];
-    }
-  }
-  A2D::Vec3 &getRefAxis() { return axis; }
-
-  void getRefAxis(TacsScalar _axis[]) {
-    _axis[0] = axis.x[0];
-    _axis[1] = axis.x[1];
-    _axis[2] = axis.x[2];
-  }
-
- private:
-  A2D::Vec3 axis;
-  /* Tolerance for colinearity test in between beam axis and ref axis */
-  const double SMALL_NUM = 1e-8;
-};
-
 /*
   The position in the beam is parametrized using the coordinates d =
   (xi, z1, z2) and the n1(xi) and n2(xi) coordinates that are
@@ -198,7 +85,7 @@ class TACSBeamRefAxisTransform : public TACSBeamTransform {
   .              u0,xi * e1^{T} * x,d0^{-1},z1 ) * T^{T} * e1
 
 
-  Note that the second coefficient for d1x can be simplifed to
+  Note that the second coefficient for d1x can be simplified to
 
   e1^{T} * x,d0^{-1},z1 * T^{T} * e1
   .   = - e1^{T} * x,d0^{-1} * x,d,z1 * x,d0^{-1} * T^{T} * e1
@@ -209,7 +96,7 @@ class TACSBeamRefAxisTransform : public TACSBeamTransform {
   d2x = T^{T} * (d2,xi * e1^{T} * x,d0^{-1} +
   .              u0,xi * e1^{T} * x,d0^{-1},z2 ) * T^{T} * e1
 
-  These expressions can be simplifed to the following:
+  These expressions can be simplified to the following:
 
   d1x = s0 * T^{T} * (d1,xi - sz1 * u0,xi)
   d2x = s0 * T^{T} * (d2,xi - sz2 * u0,xi)
@@ -223,7 +110,7 @@ class TACSBeamRefAxisTransform : public TACSBeamTransform {
 template <class quadrature, class basis, class director, class model>
 class TACSBeamElement : public TACSElement {
  public:
-  // Offset within the solution vector to the roational
+  // Offset within the solution vector to the rotational
   // parametrization defined via the director class. Here the offset
   // is 3 corresponding to the (u, v, w) displacements of the beam.
   static const int offset = 3;
@@ -310,14 +197,14 @@ class TACSBeamElement : public TACSElement {
 
   TACSElement *createElementInertialForce(const TacsScalar inertiaVec[]) {
     return new TACSBeamInertialForce<vars_per_node, quadrature, basis>(
-        con, inertiaVec);
+        transform, con, inertiaVec);
   }
 
   TACSElement *createElementCentrifugalForce(const TacsScalar omegaVec[],
                                              const TacsScalar rotCenter[],
                                              const bool first_order = false) {
     return new TACSBeamCentrifugalForce<vars_per_node, quadrature, basis>(
-        con, omegaVec, rotCenter);
+        transform, con, omegaVec, rotCenter);
   }
 
   void computeEnergies(int elemIndex, double time, const TacsScalar Xpts[],
@@ -1429,11 +1316,14 @@ int TACSBeamElement<quadrature, basis, director, model>::evalPointQuantity(
     return 3;
   } else if (quantityType == TACS_ELEMENT_DENSITY_MOMENT) {
     if (quantity) {
-      TacsScalar density = con->evalDensity(elemIndex, pt, X0.x);
+      TacsScalar mass_moment[6];
+      con->evalMassMoments(elemIndex, pt, X0.x, mass_moment);
+      TacsScalar density = mass_moment[0];
 
-      quantity[0] = density * X0.x[0];
-      quantity[1] = density * X0.x[1];
-      quantity[2] = density * X0.x[2];
+      for (int i = 0; i < 3; i++) {
+        quantity[i] = density * X0.x[i] - mass_moment[1] * n1.x[i] -
+                      mass_moment[2] * n2.x[i];
+      }
     }
 
     return 3;
@@ -1444,21 +1334,25 @@ int TACSBeamElement<quadrature, basis, director, model>::evalPointQuantity(
       // Evaluate the self MOI
       TacsScalar moments[6];
       con->evalMassMoments(elemIndex, pt, X0.x, moments);
-      I0[3] = moments[3];
-      I0[4] = moments[5];
-      I0[5] = moments[4];
+      TacsScalar density = moments[0];
+      I0[3] = moments[4] - moments[2] * moments[2] / density;
+      I0[4] = -moments[5] + moments[1] * moments[2] / density;
+      I0[5] = moments[3] - moments[1] * moments[1] / density;
       // Compute T*I0*T^{T}
       mat3x3SymmTransform(T.A, I0, quantity);
-
-      TacsScalar density = con->evalDensity(elemIndex, pt, X0.x);
+      TacsScalar dXcg[3];
+      for (int i = 0; i < 3; i++) {
+        dXcg[i] =
+            X0.x[i] - (moments[1] * n1.x[i] + moments[2] * n2.x[i]) / density;
+      }
 
       // Use parallel axis theorem to move MOI to origin
-      quantity[0] += density * (X0.x[1] * X0.x[1] + X0.x[2] * X0.x[2]);
-      quantity[1] += -density * X0.x[0] * X0.x[1];
-      quantity[2] += -density * X0.x[0] * X0.x[2];
-      quantity[3] += density * (X0.x[0] * X0.x[0] + X0.x[2] * X0.x[2]);
-      quantity[4] += -density * X0.x[2] * X0.x[1];
-      quantity[5] += density * (X0.x[0] * X0.x[0] + X0.x[1] * X0.x[1]);
+      quantity[0] += density * (dXcg[1] * dXcg[1] + dXcg[2] * dXcg[2]);
+      quantity[1] += -density * dXcg[0] * dXcg[1];
+      quantity[2] += -density * dXcg[0] * dXcg[2];
+      quantity[3] += density * (dXcg[0] * dXcg[0] + dXcg[2] * dXcg[2]);
+      quantity[4] += -density * dXcg[2] * dXcg[1];
+      quantity[5] += density * (dXcg[0] * dXcg[0] + dXcg[1] * dXcg[1]);
     }
 
     return 6;
@@ -1660,36 +1554,100 @@ void TACSBeamElement<quadrature, basis, director, model>::
     con->addStressDVSens(elemIndex, scale * dfdq[0], pt, X0.x, e, e, dvLen,
                          dfdx);
   } else if (quantityType == TACS_ELEMENT_DENSITY_MOMENT) {
-    TacsScalar dfdm = 0.0;
+    TacsScalar dfdmom[6] = {0.0};
 
     for (int i = 0; i < 3; i++) {
-      dfdm += scale * dfdq[i] * X0.x[i];
+      dfdmom[0] += scale * dfdq[i] * X0.x[i];
+      dfdmom[1] += -scale * dfdq[i] * n1.x[i];
+      dfdmom[2] += -scale * dfdq[i] * n2.x[i];
     }
 
-    con->addDensityDVSens(elemIndex, dfdm, pt, X0.x, dvLen, dfdx);
+    con->addMassMomentsDVSens(elemIndex, pt, X0.x, dfdmom, dvLen, dfdx);
   } else if (quantityType == TACS_ELEMENT_MOMENT_OF_INERTIA) {
     TacsScalar dfdI0[6] = {0.0};
 
     // Evaluate the self MOI
+    TacsScalar moments[6];
+    con->evalMassMoments(elemIndex, pt, X0.x, moments);
+    TacsScalar density = moments[0];
+
     TacsScalar dfdmoments[6] = {0.0};
     mat3x3SymmTransformSens(T.A, dfdq, dfdI0);
-    dfdmoments[3] = scale * dfdI0[3];
-    dfdmoments[5] = scale * dfdI0[4];
-    dfdmoments[4] = scale * dfdI0[5];
-
-    con->addMassMomentsDVSens(elemIndex, pt, X0.x, dfdmoments, dvLen, dfdx);
-
-    TacsScalar dfdm = 0.0;
+    dfdmoments[0] = scale *
+                    (moments[2] * moments[2] * dfdI0[3] +
+                     -moments[1] * moments[2] * dfdI0[4] +
+                     moments[1] * moments[1] * dfdI0[5]) /
+                    density / density;
+    dfdmoments[1] = scale *
+                    (-2.0 * moments[1] * dfdI0[5] + moments[2] * dfdI0[4]) /
+                    density;
+    dfdmoments[2] = scale *
+                    (-2.0 * moments[2] * dfdI0[3] + moments[1] * dfdI0[4]) /
+                    density;
+    dfdmoments[3] = scale * dfdI0[5];
+    dfdmoments[4] = scale * dfdI0[3];
+    dfdmoments[5] = -scale * dfdI0[4];
+
+    TacsScalar dXcg[3], dXcgdrho[3], dXcgdmom1[3], dXcgdmom2[3];
+    for (int i = 0; i < 3; i++) {
+      dXcg[i] =
+          X0.x[i] - (moments[1] * n1.x[i] + moments[2] * n2.x[i]) / density;
+      dXcgdrho[i] =
+          (moments[1] * n1.x[i] + moments[2] * n2.x[i]) / density / density;
+      dXcgdmom2[i] = -n2.x[i] / density;
+      dXcgdmom1[i] = -n1.x[i] / density;
+    }
 
     // Use parallel axis theorem to move MOI to origin
-    dfdm += scale * dfdq[0] * (X0.x[1] * X0.x[1] + X0.x[2] * X0.x[2]);
-    dfdm -= scale * dfdq[1] * X0.x[0] * X0.x[1];
-    dfdm -= scale * dfdq[2] * X0.x[0] * X0.x[2];
-    dfdm += scale * dfdq[3] * (X0.x[0] * X0.x[0] + X0.x[2] * X0.x[2]);
-    dfdm -= scale * dfdq[4] * X0.x[2] * X0.x[1];
-    dfdm += scale * dfdq[5] * (X0.x[0] * X0.x[0] + X0.x[1] * X0.x[1]);
-
-    con->addDensityDVSens(elemIndex, dfdm, pt, X0.x, dvLen, dfdx);
+    dfdmoments[0] +=
+        scale * dfdq[0] *
+        (dXcg[1] * dXcg[1] + dXcg[2] * dXcg[2] +
+         2.0 * density * (dXcg[1] * dXcgdrho[1] + dXcg[2] * dXcgdrho[2]));
+    dfdmoments[0] -= scale * dfdq[1] *
+                     (dXcg[0] * dXcg[1] + density * (dXcgdrho[0] * dXcg[1] +
+                                                     dXcg[0] * dXcgdrho[1]));
+    dfdmoments[0] -= scale * dfdq[2] *
+                     (dXcg[0] * dXcg[2] + density * (dXcgdrho[0] * dXcg[2] +
+                                                     dXcg[0] * dXcgdrho[2]));
+    dfdmoments[0] +=
+        scale * dfdq[3] *
+        (dXcg[0] * dXcg[0] + dXcg[2] * dXcg[2] +
+         2.0 * density * (dXcg[0] * dXcgdrho[0] + dXcg[2] * dXcgdrho[2]));
+    dfdmoments[0] -= scale * dfdq[4] *
+                     (dXcg[2] * dXcg[1] + density * (dXcgdrho[1] * dXcg[2] +
+                                                     dXcg[1] * dXcgdrho[2]));
+    dfdmoments[0] +=
+        scale * dfdq[5] *
+        (dXcg[0] * dXcg[0] + dXcg[1] * dXcg[1] +
+         2.0 * density * (dXcg[0] * dXcgdrho[0] + dXcg[1] * dXcgdrho[1]));
+
+    dfdmoments[1] += scale * dfdq[0] * density * 2.0 *
+                     (dXcg[1] * dXcgdmom1[1] + dXcg[2] * dXcgdmom1[2]);
+    dfdmoments[1] += -scale * dfdq[1] * density *
+                     (dXcgdmom1[0] * dXcg[1] + dXcg[0] * dXcgdmom1[1]);
+    dfdmoments[1] += -scale * dfdq[2] * density *
+                     (dXcgdmom1[0] * dXcg[2] + dXcg[0] * dXcgdmom1[2]);
+    dfdmoments[1] += scale * dfdq[3] * density * 2.0 *
+                     (dXcg[0] * dXcgdmom1[0] + dXcg[2] * dXcgdmom1[2]);
+    dfdmoments[1] += -scale * dfdq[4] * density *
+                     (dXcgdmom1[2] * dXcg[1] + dXcg[2] * dXcgdmom1[1]);
+    dfdmoments[1] += scale * dfdq[5] * density * 2.0 *
+                     (dXcg[0] * dXcgdmom1[0] + dXcg[1] * dXcgdmom1[1]);
+
+    dfdmoments[2] += scale * dfdq[0] * density * 2.0 *
+                     (dXcg[1] * dXcgdmom2[1] + dXcg[2] * dXcgdmom2[2]);
+    dfdmoments[2] += -scale * dfdq[1] * density *
+                     (dXcgdmom2[0] * dXcg[1] + dXcg[0] * dXcgdmom2[1]);
+    dfdmoments[2] += -scale * dfdq[2] * density *
+                     (dXcgdmom2[0] * dXcg[2] + dXcg[0] * dXcgdmom2[2]);
+    dfdmoments[2] += scale * dfdq[3] * density * 2.0 *
+                     (dXcg[0] * dXcgdmom2[0] + dXcg[2] * dXcgdmom2[2]);
+    dfdmoments[2] += -scale * dfdq[4] * density *
+                     (dXcgdmom2[2] * dXcg[1] + dXcg[2] * dXcgdmom2[1]);
+    dfdmoments[2] += scale * dfdq[5] * density * 2.0 *
+                     (dXcg[0] * dXcgdmom2[0] + dXcg[1] * dXcgdmom2[1]);
+
+    con->addMassMomentsDVSens(elemIndex, pt, X0.x, dfdmoments, dvLen, dfdx);
   }
 }
 
@@ -2023,11 +1981,15 @@ void TACSBeamElement<quadrature, basis, director, model>::
     }
   } else if (quantityType == TACS_ELEMENT_DENSITY_MOMENT) {
     // Compute the sensitivity of the strain energy density w.r.t. the strain
-    TacsScalar density = con->evalDensity(elemIndex, pt, X0.x);
+    TacsScalar mass_moment[6];
+    con->evalMassMoments(elemIndex, pt, X0.x, mass_moment);
+    TacsScalar density = mass_moment[0];
 
-    X0.xd[0] = density * dfdq[0];
-    X0.xd[1] = density * dfdq[1];
-    X0.xd[2] = density * dfdq[2];
+    for (int i = 0; i < 3; i++) {
+      X0.xd[i] = density * dfdq[i];
+      n1.xd[i] = mass_moment[2] * dfdq[i];
+      n2.xd[i] = -mass_moment[1] * dfdq[i];
+    }
   } else if (quantityType == TACS_ELEMENT_MOMENT_OF_INERTIA) {
     TACSElement::addPointQuantityXptSens(elemIndex, quantityType, time, scale,
                                          n, pt, Xpts, vars, dvars, ddvars,
diff --git a/src/elements/shell/TACSBeamElementTransform.h b/src/elements/shell/TACSBeamElementTransform.h
new file mode 100644
index 000000000..3391c5ef4
--- /dev/null
+++ b/src/elements/shell/TACSBeamElementTransform.h
@@ -0,0 +1,121 @@
+#ifndef TACS_BEAM_ELEMENT_TRANSFORM_H
+#define TACS_BEAM_ELEMENT_TRANSFORM_H
+
+#include "TACSElement.h"
+#include "a2d.h"
+
+/*
+  Compute the transformation from the local coordinates
+*/
+class TACSBeamTransform : public TACSObject {
+ public:
+  /*
+    Given the local beam element reference frame Xf, compute the
+    transformation from the global coordinates to the shell-aligned local axis.
+  */
+  virtual void computeTransform(const TacsScalar Xxi[], TacsScalar T[]) = 0;
+  virtual void addTransformSens(const TacsScalar X0xi_vals[],
+                                const TacsScalar dTvals[],
+                                TacsScalar dX0xi[]) = 0;
+  virtual A2D::Vec3 &getRefAxis() = 0;
+};
+
+/*
+  Compute the transformation
+*/
+class TACSBeamRefAxisTransform : public TACSBeamTransform {
+ public:
+  TACSBeamRefAxisTransform(const TacsScalar axis_dir[]) {
+    A2D::Vec3 axdir(axis_dir);
+    A2D::Vec3Normalize normalize(axdir, axis);
+  }
+
+  void computeTransform(const TacsScalar X0xi_vals[], TacsScalar Tvals[]) {
+    // Normalize the first direction.
+    A2D::Vec3 X0xi(X0xi_vals);
+    A2D::Vec3 t1;
+    A2D::Vec3Normalize normalizet1(X0xi, t1);
+
+    // t2_dir = axis - dot(t1, axis) * t1
+    A2D::Vec3 t2_dir;
+    A2D::Scalar dot;
+    A2D::Vec3Dot dott1(axis, t1, dot);
+    A2D::Vec3Axpy axpy(-1.0, dot, t1, axis, t2_dir);
+
+    // Check if ref axis is parallel to beam
+    if (abs(TacsRealPart(dot.value)) > 1.0 - SMALL_NUM) {
+      fprintf(stderr,
+              "TACSBeamRefAxisTransform: Error, user-provided reference axis "
+              "is parallel to beam axis. "
+              "Element behavior may be ill-conditioned.\n");
+    }
+
+    // Compute the t2 direction
+    A2D::Vec3 t2;
+    A2D::Vec3Normalize normalizet2(t2_dir, t2);
+
+    // Compute the n2 direction
+    A2D::Vec3 t3;
+    A2D::Vec3CrossProduct cross(t1, t2, t3);
+
+    // Assemble the reference frame
+    A2D::Mat3x3 T;
+    A2D::Mat3x3FromThreeVec3 assembleT(t1, t2, t3, T);
+
+    for (int i = 0; i < 9; i++) {
+      Tvals[i] = T.A[i];
+    }
+  }
+
+  void addTransformSens(const TacsScalar X0xi_vals[], const TacsScalar dTvals[],
+                        TacsScalar dX0xi[]) {
+    // Normalize the first direction.
+    A2D::ADVec3 X0xi(X0xi_vals);
+    A2D::ADVec3 t1;
+    A2D::ADVec3Normalize normalizet1(X0xi, t1);
+
+    // t2_dir = axis - dot(t1, axis) * t1
+    A2D::ADVec3 t2_dir;
+    A2D::ADScalar dot;
+    A2D::Vec3ADVecDot dott1(axis, t1, dot);
+    A2D::ADVec3VecADScalarAxpy axpy(-1.0, dot, t1, axis, t2_dir);
+
+    // Compute the t2 direction
+    A2D::ADVec3 t2;
+    A2D::ADVec3Normalize normalizet2(t2_dir, t2);
+
+    // Compute the n2 direction
+    A2D::ADVec3 t3;
+    A2D::ADVec3CrossProduct cross(t1, t2, t3);
+
+    // Assemble the referece frame
+    A2D::ADMat3x3 T(NULL, dTvals);  // Set the seeds for T
+    A2D::ADMat3x3FromThreeADVec3 assembleT(t1, t2, t3, T);
+
+    // Reverse the operations to get the derivative w.r.t. X0
+    assembleT.reverse();
+    cross.reverse();
+    normalizet2.reverse();
+    axpy.reverse();
+    dott1.reverse();
+    normalizet1.reverse();
+
+    for (int i = 0; i < 3; i++) {
+      dX0xi[i] += X0xi.xd[i];
+    }
+  }
+  A2D::Vec3 &getRefAxis() { return axis; }
+
+  void getRefAxis(TacsScalar _axis[]) {
+    _axis[0] = axis.x[0];
+    _axis[1] = axis.x[1];
+    _axis[2] = axis.x[2];
+  }
+
+ private:
+  A2D::Vec3 axis;
+  /* Tolerance for colinearity test in between beam axis and ref axis */
+  const double SMALL_NUM = 1e-8;
+};
+
+#endif  // TACS_BEAM_ELEMENT_TRANSFORM_H
\ No newline at end of file
diff --git a/src/elements/shell/TACSBeamInertialForce.h b/src/elements/shell/TACSBeamInertialForce.h
index 75f6b21ee..913d03ce3 100644
--- a/src/elements/shell/TACSBeamInertialForce.h
+++ b/src/elements/shell/TACSBeamInertialForce.h
@@ -4,6 +4,7 @@
 #include "TACSBeamConstitutive.h"
 #include "TACSBeamElementBasis.h"
 #include "TACSBeamElementQuadrature.h"
+#include "TACSBeamElementTransform.h"
 #include "TACSBeamUtilities.h"
 #include "TACSElement.h"
 #include "TACSElementAlgebra.h"
@@ -14,14 +15,22 @@
 template <int vars_per_node, class quadrature, class basis>
 class TACSBeamInertialForce : public TACSElement {
  public:
-  TACSBeamInertialForce(TACSBeamConstitutive *_con,
+  TACSBeamInertialForce(TACSBeamTransform *_transform,
+                        TACSBeamConstitutive *_con,
                         const TacsScalar _inertiaVec[]) {
+    transform = _transform;
+    transform->incref();
+
     con = _con;
     con->incref();
+
     memcpy(inertiaVec, _inertiaVec, 3 * sizeof(TacsScalar));
   }
 
   ~TACSBeamInertialForce() {
+    if (transform) {
+      transform->decref();
+    }
     if (con) {
       con->decref();
     }
@@ -78,6 +87,13 @@ class TACSBeamInertialForce : public TACSElement {
     // Compute the number of quadrature points
     const int nquad = quadrature::getNumQuadraturePoints();
 
+    // Get the reference axis
+    const A2D::Vec3 &axis = transform->getRefAxis();
+
+    // Compute the normal directions
+    TacsScalar fn1[3 * basis::NUM_NODES], fn2[3 * basis::NUM_NODES];
+    TacsBeamComputeNodeNormals<basis>(Xpts, axis, fn1, fn2);
+
     // Loop over each quadrature point and add the residual contribution
     for (int quad_index = 0; quad_index < nquad; quad_index++) {
       // Get the quadrature weight
@@ -87,28 +103,42 @@ class TACSBeamInertialForce : public TACSElement {
       // Tangent to the beam
       A2D::Vec3 X0, X0xi;
 
+      // Interpolated normal directions
+      A2D::ADVec3 n1, n2;
+
       // Compute X, X,xi and the interpolated normal
       basis::template interpFields<3, 3>(pt, Xpts, X0.x);
       basis::template interpFieldsGrad<3, 3>(pt, Xpts, X0xi.x);
+      basis::template interpFields<3, 3>(pt, fn1, n1.x);
+      basis::template interpFields<3, 3>(pt, fn2, n2.x);
 
       // Compute the determinant of the transform
       A2D::Scalar detXd;
       A2D::Vec3Norm(X0xi, detXd);
 
-      TacsScalar mass = con->evalDensity(elemIndex, pt, X0.x);
+      TacsScalar mass_moment[6];
+      con->evalMassMoments(elemIndex, pt, X0.x, mass_moment);
+      TacsScalar mass = mass_moment[0];
 
       // Compute the traction
-      TacsScalar tr[3];
+      TacsScalar tr[vars_per_node] = {0.0};
       tr[0] = -detXd.value * weight * mass * inertiaVec[0];
       tr[1] = -detXd.value * weight * mass * inertiaVec[1];
       tr[2] = -detXd.value * weight * mass * inertiaVec[2];
-
-      basis::template addInterpFieldsTranspose<vars_per_node, 3>(pt, tr, res);
+      // Add moment terms if theres a beam offset
+      crossProductAdd(detXd.value * weight * mass_moment[1], n1.x, inertiaVec,
+                      &tr[3]);
+      crossProductAdd(detXd.value * weight * mass_moment[2], n2.x, inertiaVec,
+                      &tr[3]);
+
+      basis::template addInterpFieldsTranspose<vars_per_node, vars_per_node>(
+          pt, tr, res);
     }
   }
 
  private:
   TacsScalar inertiaVec[3];
+  TACSBeamTransform *transform;
   TACSBeamConstitutive *con;
 };
 
diff --git a/src/elements/shell/TACSShellCentrifugalForce.h b/src/elements/shell/TACSShellCentrifugalForce.h
index 3790b0230..a726fd31a 100644
--- a/src/elements/shell/TACSShellCentrifugalForce.h
+++ b/src/elements/shell/TACSShellCentrifugalForce.h
@@ -104,7 +104,9 @@ class TACSShellCentrifugalForce : public TACSElement {
       TacsScalar detXd = det3x3(Xd);
       detXd *= weight;
 
-      TacsScalar mass = con->evalDensity(elemIndex, pt, X);
+      TacsScalar moments[3];
+      con->evalMassMoments(elemIndex, pt, X, moments);
+      TacsScalar mass = moments[0];
 
       TacsScalar r[3], wxr[3], ac[3];
 
@@ -113,6 +115,13 @@ class TACSShellCentrifugalForce : public TACSElement {
       r[1] = X[1] - rotCenter[1] + U[1];
       r[2] = X[2] - rotCenter[2] + U[2];
 
+      // Account for shell mass offset
+      if (!first_order) {
+        for (int i = 0; i < 3; i++) {
+          r[i] -= moments[1] / mass * n[i];
+        }
+      }
+
       // Compute omega x r
       crossProduct(omegaVec, r, wxr);
 
@@ -120,12 +129,17 @@ class TACSShellCentrifugalForce : public TACSElement {
       crossProduct(omegaVec, wxr, ac);
 
       // Compute the traction
-      TacsScalar tr[3];
+      TacsScalar tr[vars_per_node] = {0.0};
       tr[0] = detXd * mass * ac[0];
       tr[1] = detXd * mass * ac[1];
       tr[2] = detXd * mass * ac[2];
+      // Add moment terms if theres a shell offset
+      if (!first_order) {
+        crossProductAdd(-detXd * moments[1], n, ac, &tr[3]);
+      }
 
-      basis::template addInterpFieldsTranspose<vars_per_node, 3>(pt, tr, res);
+      basis::template addInterpFieldsTranspose<vars_per_node, vars_per_node>(
+          pt, tr, res);
     }
   }
 
diff --git a/src/elements/shell/TACSShellElement.h b/src/elements/shell/TACSShellElement.h
index aafe630a0..8687348b5 100644
--- a/src/elements/shell/TACSShellElement.h
+++ b/src/elements/shell/TACSShellElement.h
@@ -916,11 +916,13 @@ int TACSShellElement<quadrature, basis, director, model>::evalPointQuantity(
       TacsScalar Xd[9];
       TacsShellAssembleFrame(Xxi, n0, Xd);
       *detXd = det3x3(Xd);
-      TacsScalar density = con->evalDensity(elemIndex, pt, X);
+      TacsScalar moments[3];
+      con->evalMassMoments(elemIndex, pt, X, moments);
+      TacsScalar density = moments[0];
 
-      quantity[0] = density * X[0];
-      quantity[1] = density * X[1];
-      quantity[2] = density * X[2];
+      quantity[0] = density * X[0] + moments[1] * n0[0];
+      quantity[1] = density * X[1] + moments[1] * n0[1];
+      quantity[2] = density * X[2] + moments[1] * n0[2];
     }
 
     return 3;
@@ -938,24 +940,28 @@ int TACSShellElement<quadrature, basis, director, model>::evalPointQuantity(
       TacsScalar Xd[9];
       TacsShellAssembleFrame(Xxi, n0, Xd);
       *detXd = det3x3(Xd);
-      TacsScalar density = con->evalDensity(elemIndex, pt, X);
 
       TacsScalar I0[6] = {0.0};
 
       // Evaluate the self MOI
       TacsScalar moments[3];
       con->evalMassMoments(elemIndex, pt, X, moments);
-      I0[0] = I0[3] = moments[2];
+      TacsScalar density = moments[0];
+      I0[0] = I0[3] = moments[2] - moments[1] * moments[1] / density;
       // Compute T*I0*T^{T}
       mat3x3SymmTransform(T, I0, quantity);
+      TacsScalar dXcg[3];
+      for (int i = 0; i < 3; i++) {
+        dXcg[i] = X[i] + moments[1] / density * n0[i];
+      }
 
       // Use parallel axis theorem to move MOI to origin
-      quantity[0] += density * (X[1] * X[1] + X[2] * X[2]);
-      quantity[1] += -density * X[0] * X[1];
-      quantity[2] += -density * X[0] * X[2];
-      quantity[3] += density * (X[0] * X[0] + X[2] * X[2]);
-      quantity[4] += -density * X[2] * X[1];
-      quantity[5] += density * (X[0] * X[0] + X[1] * X[1]);
+      quantity[0] += density * (dXcg[1] * dXcg[1] + dXcg[2] * dXcg[2]);
+      quantity[1] += -density * dXcg[0] * dXcg[1];
+      quantity[2] += -density * dXcg[0] * dXcg[2];
+      quantity[3] += density * (dXcg[0] * dXcg[0] + dXcg[2] * dXcg[2]);
+      quantity[4] += -density * dXcg[2] * dXcg[1];
+      quantity[5] += density * (dXcg[0] * dXcg[0] + dXcg[1] * dXcg[1]);
     }
 
     return 6;
@@ -1050,16 +1056,23 @@ void TACSShellElement<quadrature, basis, director, model>::
 
     con->addDensityDVSens(elemIndex, scale * dfdq[0], pt, X, dvLen, dfdx);
   } else if (quantityType == TACS_ELEMENT_DENSITY_MOMENT) {
-    TacsScalar X[3];
+    // Compute the node normal directions
+    TacsScalar fn[3 * num_nodes];
+    TacsShellComputeNodeNormals<basis>(Xpts, fn);
+
+    // Compute X and the interpolated normal n0
+    TacsScalar X[3], n0[3];
     basis::template interpFields<3, 3>(pt, Xpts, X);
+    basis::template interpFields<3, 3>(pt, fn, n0);
 
-    TacsScalar dfdm = 0.0;
+    TacsScalar dfdmoments[3] = {0.0};
 
     for (int i = 0; i < 3; i++) {
-      dfdm += scale * dfdq[i] * X[i];
+      dfdmoments[0] += scale * dfdq[i] * X[i];
+      dfdmoments[1] += scale * dfdq[i] * n0[i];
     }
 
-    con->addDensityDVSens(elemIndex, dfdm, pt, X, dvLen, dfdx);
+    con->addMassMomentsDVSens(elemIndex, pt, X, dfdmoments, dvLen, dfdx);
   } else if (quantityType == TACS_ELEMENT_MOMENT_OF_INERTIA) {
     // Compute the node normal directions
     TacsScalar fn[3 * num_nodes];
@@ -1080,23 +1093,60 @@ void TACSShellElement<quadrature, basis, director, model>::
     TacsScalar dfdI0[6] = {0.0};
 
     // Evaluate the self MOI
-    TacsScalar dfdmoments[3] = {0.0};
+    TacsScalar moments[3];
+    con->evalMassMoments(elemIndex, pt, X, moments);
+    TacsScalar density = moments[0];
+
+    // Evaluate the self MOI
+    TacsScalar dfdmoments[3];
     mat3x3SymmTransformSens(T, dfdq, dfdI0);
     dfdmoments[2] = scale * (dfdI0[0] + dfdI0[3]);
+    dfdmoments[1] = -scale * 2.0 * moments[1] / density * (dfdI0[0] + dfdI0[3]);
+    dfdmoments[0] = scale * moments[1] * moments[1] / density / density *
+                    (dfdI0[0] + dfdI0[3]);
 
-    con->addMassMomentsDVSens(elemIndex, pt, X, dfdmoments, dvLen, dfdx);
-
-    TacsScalar dfdm = 0.0;
+    TacsScalar dXcg[3];
+    for (int i = 0; i < 3; i++) {
+      dXcg[i] = X[i] + moments[1] / density * n0[i];
+    }
 
     // Use parallel axis theorem to move MOI to origin
-    dfdm += scale * dfdq[0] * (X[1] * X[1] + X[2] * X[2]);
-    dfdm -= scale * dfdq[1] * X[0] * X[1];
-    dfdm -= scale * dfdq[2] * X[0] * X[2];
-    dfdm += scale * dfdq[3] * (X[0] * X[0] + X[2] * X[2]);
-    dfdm -= scale * dfdq[4] * X[2] * X[1];
-    dfdm += scale * dfdq[5] * (X[0] * X[0] + X[1] * X[1]);
-
-    con->addDensityDVSens(elemIndex, dfdm, pt, X, dvLen, dfdx);
+    dfdmoments[0] +=
+        scale * dfdq[0] *
+        (dXcg[1] * dXcg[1] + dXcg[2] * dXcg[2] -
+         2.0 * moments[1] / density * (dXcg[1] * n0[1] + dXcg[2] * n0[2]));
+    dfdmoments[0] -=
+        scale * dfdq[1] *
+        (dXcg[0] * dXcg[1] -
+         moments[1] / density * (dXcg[0] * n0[1] + dXcg[1] * n0[0]));
+    dfdmoments[0] -=
+        scale * dfdq[2] *
+        (dXcg[0] * dXcg[2] -
+         moments[1] / density * (dXcg[0] * n0[2] + dXcg[2] * n0[0]));
+    dfdmoments[0] +=
+        scale * dfdq[3] *
+        (dXcg[0] * dXcg[0] + dXcg[2] * dXcg[2] -
+         2.0 * moments[1] / density * (dXcg[0] * n0[0] + dXcg[2] * n0[2]));
+    dfdmoments[0] -=
+        scale * dfdq[4] *
+        (dXcg[2] * dXcg[1] -
+         moments[1] / density * (dXcg[1] * n0[2] + dXcg[2] * n0[1]));
+    dfdmoments[0] +=
+        scale * dfdq[5] *
+        (dXcg[0] * dXcg[0] + dXcg[1] * dXcg[1] -
+         2.0 * moments[1] / density * (dXcg[0] * n0[0] + dXcg[1] * n0[1]));
+
+    dfdmoments[1] +=
+        scale * dfdq[0] * 2.0 * (dXcg[1] * n0[1] + dXcg[2] * n0[2]);
+    dfdmoments[1] -= scale * dfdq[1] * (dXcg[0] * n0[1] + dXcg[1] * n0[0]);
+    dfdmoments[1] -= scale * dfdq[2] * (dXcg[0] * n0[2] + dXcg[2] * n0[0]);
+    dfdmoments[1] +=
+        scale * dfdq[3] * 2.0 * (dXcg[0] * n0[0] + dXcg[2] * n0[2]);
+    dfdmoments[1] -= scale * dfdq[4] * (dXcg[1] * n0[2] + dXcg[2] * n0[1]);
+    dfdmoments[1] +=
+        scale * dfdq[5] * 2.0 * (dXcg[0] * n0[0] + dXcg[1] * n0[1]);
+
+    con->addMassMomentsDVSens(elemIndex, pt, X, dfdmoments, dvLen, dfdx);
   }
 }
 
diff --git a/src/elements/shell/TACSShellInertialForce.h b/src/elements/shell/TACSShellInertialForce.h
index 33f4c1996..c8372fd64 100644
--- a/src/elements/shell/TACSShellInertialForce.h
+++ b/src/elements/shell/TACSShellInertialForce.h
@@ -96,15 +96,20 @@ class TACSShellInertialForce : public TACSElement {
       TacsScalar detXd = det3x3(Xd);
       detXd *= weight;
 
-      TacsScalar mass = con->evalDensity(elemIndex, pt, X);
+      TacsScalar moments[3];
+      con->evalMassMoments(elemIndex, pt, X, moments);
+      TacsScalar mass = moments[0];
 
       // Compute the traction
-      TacsScalar tr[3];
+      TacsScalar tr[vars_per_node] = {0.0};
       tr[0] = -detXd * mass * inertiaVec[0];
       tr[1] = -detXd * mass * inertiaVec[1];
       tr[2] = -detXd * mass * inertiaVec[2];
+      // Add moment terms if theres a shell offset
+      crossProductAdd(detXd * moments[1], n, inertiaVec, &tr[3]);
 
-      basis::template addInterpFieldsTranspose<vars_per_node, 3>(pt, tr, res);
+      basis::template addInterpFieldsTranspose<vars_per_node, vars_per_node>(
+          pt, tr, res);
     }
   }
 
@@ -113,4 +118,4 @@ class TACSShellInertialForce : public TACSElement {
   TACSShellConstitutive *con;
 };
 
-#endif  // TACS_SHELL_TRACTION_H
\ No newline at end of file
+#endif  // TACS_SHELL_INERTIAL_FORCE_H
\ No newline at end of file
diff --git a/tacs/constitutive.pyx b/tacs/constitutive.pyx
index 96c810d5d..1d97bd1d3 100644
--- a/tacs/constitutive.pyx
+++ b/tacs/constitutive.pyx
@@ -671,6 +671,9 @@ cdef class IsoShellConstitutive(ShellConstitutive):
             (i.e. no design variable).
         tlb (float or complex, optional): Thickness variable lower bound (keyword argument). Defaults to 0.0.
         tub (float or complex, optional): Thickness variable upper bound (keyword argument). Defaults to 10.0.
+        tOffset (float or complex, optional): Offset distance of reference plane (where nodes are located) relative to thickness mid-plane.
+            Measured in fraction of shell thickness. A value of 0.5 places the reference plane at the top of the plate,
+            a value of 0.0 at the plate mid-plane, and a value of -0.5 at the bottom of the plate. Defaults to 0.0.
     """
     def __cinit__(self, *args, **kwargs):
         cdef TACSMaterialProperties *props = NULL
@@ -678,6 +681,7 @@ cdef class IsoShellConstitutive(ShellConstitutive):
         cdef int tNum = -1
         cdef TacsScalar tlb = 0.0
         cdef TacsScalar tub = 10.0
+        cdef TacsScalar tOff = 0.0
 
         if len(args) >= 1:
             props = (<MaterialProperties>args[0]).ptr
@@ -690,10 +694,12 @@ cdef class IsoShellConstitutive(ShellConstitutive):
             tlb = kwargs['tlb']
         if 'tub' in kwargs:
             tub = kwargs['tub']
+        if 'tOffset' in kwargs:
+            tOff = kwargs['tOffset']
 
         if props is not NULL:
             self.cptr = new TACSIsoShellConstitutive(props, t, tNum,
-                                                     tlb, tub)
+                                                     tlb, tub, tOff)
             self.ptr = self.cptr
             self.ptr.incref()
         else:
@@ -731,11 +737,15 @@ cdef class CompositeShellConstitutive(ShellConstitutive):
        ply_thicknesses (numpy.ndarray[float or complex]): Array of ply thicknesses in layup.
        ply_angles (numpy.ndarray[float or complex]): Array of ply angles (in radians) in layup.
        kcorr (float or complex, optional): FSDT shear correction factor. Defaults to 5.0/6.0.
+       tOffset (float or complex, optional): Offset distance of reference plane (where nodes are located) relative to thickness mid-plane.
+            Measured in fraction of shell thickness. A value of 0.5 places the reference plane at the top of the plate,
+            a value of 0.0 at the plate mid-plane, and a value of -0.5 at the bottom of the plate. Defaults to 0.0.
     """
     def __cinit__(self, ply_list,
                   np.ndarray[TacsScalar, ndim=1, mode='c'] ply_thicknesses,
                   np.ndarray[TacsScalar, ndim=1, mode='c'] ply_angles,
-                  TacsScalar kcorr=5.0/6.0):
+                  TacsScalar kcorr=5.0/6.0,
+                  TacsScalar tOffset=0.0):
 
         num_plies = len(ply_list)
         if len(ply_thicknesses) != num_plies:
@@ -754,7 +764,7 @@ cdef class CompositeShellConstitutive(ShellConstitutive):
 
         self.cptr = new TACSCompositeShellConstitutive(num_plies, plys,
                                                        <TacsScalar*>ply_thicknesses.data,
-                                                       <TacsScalar*>ply_angles.data, kcorr)
+                                                       <TacsScalar*>ply_angles.data, kcorr, tOffset)
         self.ptr = self.cptr
         self.ptr.incref()
 
@@ -783,9 +793,14 @@ cdef class CompositeShellConstitutive(ShellConstitutive):
             ply_id = self.props[i].getNastranID()
             mat_ids.append(ply_id)
 
+        # Calculate distance from ref plane to bottom ply
+        lam_thickness = sum(ply_thicknesses)
+        z0 = -(comp_ptr.getThicknessOffset() + 0.5) * lam_thickness
+
         prop = nastran_cards.properties.shell.PCOMP(self.nastranID, mat_ids,
                                                     ply_thicknesses.astype(float),
-                                                    np.rad2deg(np.real(ply_angles)))
+                                                    np.rad2deg(np.real(ply_angles)),
+                                                    z0=np.real(z0))
         return prop
 
 cdef class BladeStiffenedShellConstitutive(ShellConstitutive):
@@ -1078,6 +1093,9 @@ cdef class SmearedCompositeShellConstitutive(ShellConstitutive):
           Defaults to 0.0.
        ply_fraction_ub (numpy.ndarray[float or complex], optional): Upper bound for ply fraction design variables (keyword argument).
           Defaults to 1.0.
+       t_offset (float or complex, optional): Offset distance of reference plane (where nodes are located) relative to thickness mid-plane.
+            Measured in fraction of shell thickness. A value of 0.5 places the reference plane at the top of the plate,
+            a value of 0.0 at the plate mid-plane, and a value of -0.5 at the bottom of the plate. Defaults to 0.0.
     """
     def __cinit__(self, ply_list,
                   TacsScalar thickness,
@@ -1088,7 +1106,8 @@ cdef class SmearedCompositeShellConstitutive(ShellConstitutive):
                   TacsScalar thickness_lb=0.0,
                   TacsScalar thickness_ub=1e20,
                   np.ndarray[TacsScalar, ndim=1, mode='c'] ply_fraction_lb=None,
-                  np.ndarray[TacsScalar, ndim=1, mode='c'] ply_fraction_ub=None,):
+                  np.ndarray[TacsScalar, ndim=1, mode='c'] ply_fraction_ub=None,
+                  TacsScalar t_offset=0.0):
 
         num_plies = len(ply_list)
 
@@ -1126,7 +1145,8 @@ cdef class SmearedCompositeShellConstitutive(ShellConstitutive):
                                                               <int*>_ply_fraction_dv_nums.data,
                                                               thickness_lb, thickness_ub,
                                                               <TacsScalar *> _ply_fraction_lb.data,
-                                                              <TacsScalar*>_ply_fraction_ub.data)
+                                                              <TacsScalar*>_ply_fraction_ub.data,
+                                                              t_offset)
         self.ptr = self.cptr
         self.ptr.incref()
 
@@ -1157,10 +1177,13 @@ cdef class SmearedCompositeShellConstitutive(ShellConstitutive):
             ply_id = self.props[i].getNastranID()
             mat_ids.append(ply_id)
 
+        z0 = -(comp_ptr.getThicknessOffset() + 0.5) * t_tot
+
         prop = nastran_cards.properties.shell.PCOMP(self.nastranID, mat_ids,
                                                     ply_thicknesses.astype(float),
                                                     np.rad2deg(np.real(ply_angles)),
-                                                    lam="SMEAR")
+                                                    lam="SMEAR",
+                                                    z0=np.real(z0))
         return prop
 
 cdef class LamParamShellConstitutive(ShellConstitutive):
@@ -1317,7 +1340,7 @@ cdef class BasicBeamConstitutive(BeamConstitutive):
         A = np.real(s[0])
         J = np.real(s[1])
         Iz = np.real(s[2])
-        Iyz = np.real(s[3])
+        Iyz = -np.real(s[3])
         e[:] = 0.0
         e[3] = 1.0 / E
         e[4:5] = 1.0 / (G * A)
@@ -1325,7 +1348,8 @@ cdef class BasicBeamConstitutive(BeamConstitutive):
         Iy = np.real(s[3])
         ky = np.real(s[4])
         kz = np.real(s[5])
-        prop = nastran_cards.properties.bars.PBAR(self.nastranID, mat_id, A, Iz, Iy, Iyz, J, k1=ky, k2=kz)
+        # Nastran uses negative convention for POI's
+        prop = nastran_cards.properties.bars.PBAR(self.nastranID, mat_id, A, Iz, Iy, -Iyz, J, k1=ky, k2=kz)
         return prop
 
 cdef class IsoTubeBeamConstitutive(BeamConstitutive):
@@ -1426,6 +1450,12 @@ cdef class IsoRectangleBeamConstitutive(BeamConstitutive):
             (i.e. no design variable).
         tlb (float or complex, optional): Lower bound on thickness (keyword argument). Defaults to 0.0.
         tub (float or complex, optional): Upper bound on thickness (keyword argument). Defaults to 10.0.
+        wOffset (float or complex, optional): Offset distance along width axis of reference axis (where nodes are located) relative to elastic axis.
+            Measured in fraction of cross-section width. A value of 0.5 places the reference axis in the plane z=/2,
+            a value of 0.0 at z=0.0, and a value of -0.5 at z=-w/2 (where z is the local beam axis). Defaults to 0.0.
+        tOffset (float or complex, optional): Offset distance along thickness axis of reference axis (where nodes are located) relative to elastic axis.
+            Measured in fraction of cross-section thickness. A value of 0.5 places the reference axis in the plane y=t/2,
+            a value of 0.0 at y=0.0, and a value of -0.5 at y=-t/2 (where y is the local beam axis). Defaults to 0.0.
     """
     def __cinit__(self, *args, **kwargs):
         cdef TACSMaterialProperties *props = NULL
@@ -1437,6 +1467,8 @@ cdef class IsoRectangleBeamConstitutive(BeamConstitutive):
         cdef int tNum = -1
         cdef TacsScalar tlb = 0.0
         cdef TacsScalar tub = 10.0
+        cdef TacsScalar wOffset = 0.0
+        cdef TacsScalar tOffset = 0.0
 
         if len(args) >= 1:
             props = (<MaterialProperties>args[0]).ptr
@@ -1457,10 +1489,14 @@ cdef class IsoRectangleBeamConstitutive(BeamConstitutive):
             tlb = kwargs['tlb']
         if 'tub' in kwargs:
             tub = kwargs['tub']
+        if 'wOffset' in kwargs:
+            wOffset = kwargs['wOffset']
+        if 'tOffset' in kwargs:
+            tOffset = kwargs['tOffset']
 
         if props is not NULL:
             self.cptr = new TACSIsoRectangleBeamConstitutive(props, w, t, wNum, tNum,
-                                                             wlb, wub, tlb, tub)
+                                                             wlb, wub, tlb, tub, wOffset, tOffset)
             self.ptr = self.cptr
             self.ptr.incref()
         else:
diff --git a/tacs/cpp_headers/constitutive.pxd b/tacs/cpp_headers/constitutive.pxd
index 5aba076f5..b449a2f30 100644
--- a/tacs/cpp_headers/constitutive.pxd
+++ b/tacs/cpp_headers/constitutive.pxd
@@ -79,14 +79,15 @@ cdef extern from "TACSShellConstitutive.h":
 cdef extern from "TACSIsoShellConstitutive.h":
     cdef cppclass TACSIsoShellConstitutive(TACSShellConstitutive):
         TACSIsoShellConstitutive(TACSMaterialProperties*, TacsScalar, int,
-                                 TacsScalar, TacsScalar)
+                                 TacsScalar, TacsScalar, TacsScalar)
 
 cdef extern from "TACSCompositeShellConstitutive.h":
     cdef cppclass TACSCompositeShellConstitutive(TACSShellConstitutive):
         TACSCompositeShellConstitutive(int, TACSOrthotropicPly**, const TacsScalar*,
-                                       const TacsScalar*, TacsScalar)
+                                       const TacsScalar*, TacsScalar, TacsScalar)
         void getPlyThicknesses(TacsScalar*);
         void getPlyAngles(TacsScalar*);
+        TacsScalar getThicknessOffset();
 
 cdef extern from "TACSLamParamShellConstitutive.h":
     cdef cppclass TACSLamParamShellConstitutive(TACSShellConstitutive):
@@ -113,10 +114,12 @@ cdef extern from "TACSSmearedCompositeShellConstitutive.h":
                                        const TacsScalar*, const TacsScalar*,
                                        int, const int*,
                                        TacsScalar, TacsScalar,
-                                       const TacsScalar*, const TacsScalar*)
+                                       const TacsScalar*, const TacsScalar*,
+                                       TacsScalar)
         TacsScalar getLaminateThickness();
         void getPlyAngles(TacsScalar*);
-        void getPlyFractions(TacsScalar *);
+        void getPlyFractions(TacsScalar*);
+        TacsScalar getThicknessOffset();
 
 cdef extern from "TACSLamParamShellConstitutive.h":
     cdef cppclass TACSLamParamShellConstitutive(TACSShellConstitutive):
@@ -214,7 +217,8 @@ cdef extern from "TACSIsoTubeBeamConstitutive.h":
 cdef extern from "TACSIsoRectangleBeamConstitutive.h":
     cdef cppclass TACSIsoRectangleBeamConstitutive(TACSBeamConstitutive):
         TACSIsoRectangleBeamConstitutive(TACSMaterialProperties*, TacsScalar, TacsScalar,
-                                         int, int, TacsScalar, TacsScalar, TacsScalar, TacsScalar)
+                                         int, int, TacsScalar, TacsScalar, TacsScalar, TacsScalar,
+                                         TacsScalar, TacsScalar)
 
 cdef extern from "TACSGeneralMassConstitutive.h":
     cdef cppclass TACSGeneralMassConstitutive(TACSConstitutive):
diff --git a/tacs/cpp_headers/elements.pxd b/tacs/cpp_headers/elements.pxd
index 62f0b49ba..efa2a8159 100644
--- a/tacs/cpp_headers/elements.pxd
+++ b/tacs/cpp_headers/elements.pxd
@@ -192,7 +192,7 @@ cdef extern from "TACSShellElementTransform.h":
         TACSShellRefAxisTransform(const TacsScalar*)
         void getRefAxis(TacsScalar*)
 
-cdef extern from "TACSBeamElement.h":
+cdef extern from "TACSBeamElementTransform.h":
     cdef cppclass TACSBeamTransform(TACSObject):
         pass
 
diff --git a/tacs/pytacs.py b/tacs/pytacs.py
index c014ad43c..bf5c4dc48 100755
--- a/tacs/pytacs.py
+++ b/tacs/pytacs.py
@@ -992,17 +992,21 @@ def elemCallBack(
                 plyThicknesses = np.array(plyThicknesses, dtype=self.dtype)
                 plyAngles = np.array(plyAngles, dtype=self.dtype)
 
+                # Get the total laminate thickness
+                lamThickness = propInfo.Thickness()
+                # Get the offset distance from the ref plane to the midplane
+                tOffset = -(propInfo.z0 / lamThickness + 0.5)
+
                 if propInfo.lam is None or propInfo.lam in ["SYM", "MEM"]:
                     # Discrete laminate class (not for optimization)
                     con = tacs.constitutive.CompositeShellConstitutive(
-                        plyMats, plyThicknesses, plyAngles
+                        plyMats, plyThicknesses, plyAngles, tOffset=tOffset
                     )
 
                 elif propInfo.lam == "SMEAR":
-                    lamThickness = sum(plyThicknesses)
                     plyFractions = plyThicknesses / lamThickness
                     con = tacs.constitutive.SmearedCompositeShellConstitutive(
-                        plyMats, lamThickness, plyAngles, plyFractions
+                        plyMats, lamThickness, plyAngles, plyFractions, t_offset=tOffset
                     )
 
                 # Need to add functionality to consider only membrane in TACS for type = MEM
@@ -1042,7 +1046,8 @@ def elemCallBack(
                 area = propInfo.A
                 I1 = propInfo.i1
                 I2 = propInfo.i2
-                I12 = propInfo.i12
+                # Nastran uses negative convention for POI's
+                I12 = -propInfo.i12
                 J = propInfo.j
                 k1 = propInfo.k1
                 k2 = propInfo.k2
diff --git a/tests/constitutive_tests/test_composite_shell_constitutive.py b/tests/constitutive_tests/test_composite_shell_constitutive.py
index 3eb73fc4e..efc7ea371 100644
--- a/tests/constitutive_tests/test_composite_shell_constitutive.py
+++ b/tests/constitutive_tests/test_composite_shell_constitutive.py
@@ -91,6 +91,8 @@ def setUp(self):
         self.layup_list = [iso_layup, ortho_layup]
         self.ply_thicknesses = np.array([ply_thickness] * nplies, dtype=self.dtype)
         self.ply_angles = np.array([0.0, -45.0, 90.0], dtype=self.dtype) * DEG2RAD
+        # Distance between mid-plane and reference plane (nodes) of shell
+        self.tOffset = 0.5 # This places the nodes at the top ply
 
         # Seed random number generator in tacs for consistent test results
         elements.SeedRandomGenerator(0)
@@ -100,7 +102,7 @@ def test_constitutive_density(self):
         for layup in self.layup_list:
             with self.subTest(layup=layup):
                 con = constitutive.CompositeShellConstitutive(
-                    layup, self.ply_thicknesses, self.ply_angles
+                    layup, self.ply_thicknesses, self.ply_angles, tOffset=self.tOffset,
                 )
                 fail = constitutive.TestConstitutiveDensity(
                     con,
@@ -120,7 +122,7 @@ def test_constitutive_specific_heat(self):
         for layup in self.layup_list:
             with self.subTest(layup=layup):
                 con = constitutive.CompositeShellConstitutive(
-                    layup, self.ply_thicknesses, self.ply_angles
+                    layup, self.ply_thicknesses, self.ply_angles, tOffset=self.tOffset,
                 )
                 fail = constitutive.TestConstitutiveSpecificHeat(
                     con,
@@ -140,7 +142,7 @@ def test_constitutive_heat_flux(self):
         for layup in self.layup_list:
             with self.subTest(layup=layup):
                 con = constitutive.CompositeShellConstitutive(
-                    layup, self.ply_thicknesses, self.ply_angles
+                    layup, self.ply_thicknesses, self.ply_angles, tOffset=self.tOffset,
                 )
                 fail = constitutive.TestConstitutiveHeatFlux(
                     con,
@@ -160,7 +162,7 @@ def test_constitutive_stress(self):
         for layup in self.layup_list:
             with self.subTest(layup=layup):
                 con = constitutive.CompositeShellConstitutive(
-                    layup, self.ply_thicknesses, self.ply_angles
+                    layup, self.ply_thicknesses, self.ply_angles, tOffset=self.tOffset,
                 )
                 fail = constitutive.TestConstitutiveStress(
                     con,
@@ -180,7 +182,7 @@ def test_constitutive_thermal_strain(self):
         for layup in self.layup_list:
             with self.subTest(layup=layup):
                 con = constitutive.CompositeShellConstitutive(
-                    layup, self.ply_thicknesses, self.ply_angles
+                    layup, self.ply_thicknesses, self.ply_angles, tOffset=self.tOffset,
                 )
                 fail = constitutive.TestConstitutiveThermalStrain(
                     con,
@@ -200,7 +202,7 @@ def test_constitutive_failure(self):
         for layup in self.layup_list:
             with self.subTest(layup=layup):
                 con = constitutive.CompositeShellConstitutive(
-                    layup, self.ply_thicknesses, self.ply_angles
+                    layup, self.ply_thicknesses, self.ply_angles, tOffset=self.tOffset,
                 )
                 fail = constitutive.TestConstitutiveFailure(
                     con,
@@ -219,7 +221,7 @@ def test_constitutive_failure_strain_sens(self):
         for layup in self.layup_list:
             with self.subTest(layup=layup):
                 con = constitutive.CompositeShellConstitutive(
-                    layup, self.ply_thicknesses, self.ply_angles
+                    layup, self.ply_thicknesses, self.ply_angles, tOffset=self.tOffset,
                 )
                 fail = constitutive.TestConstitutiveFailureStrainSens(
                     con,
diff --git a/tests/constitutive_tests/test_iso_rectangle_beam_constitutive.py b/tests/constitutive_tests/test_iso_rectangle_beam_constitutive.py
index b744641c5..46ee4cc4a 100644
--- a/tests/constitutive_tests/test_iso_rectangle_beam_constitutive.py
+++ b/tests/constitutive_tests/test_iso_rectangle_beam_constitutive.py
@@ -52,7 +52,7 @@ def setUp(self):
 
         # Create stiffness (need class)
         self.con = constitutive.IsoRectangleBeamConstitutive(
-            self.props, t=t, tNum=tNum, w=w, wNum=wNum
+            self.props, t=t, tNum=tNum, w=w, wNum=wNum, wOffset=0.5, tOffset=-0.5
         )
 
         # Seed random number generator in tacs for consistent test results
diff --git a/tests/constitutive_tests/test_iso_shell_constitutive.py b/tests/constitutive_tests/test_iso_shell_constitutive.py
index 1ba62ca03..ce30614e6 100644
--- a/tests/constitutive_tests/test_iso_shell_constitutive.py
+++ b/tests/constitutive_tests/test_iso_shell_constitutive.py
@@ -13,7 +13,7 @@ def setUp(self):
             self.rtol = 1e-11
         else:
             self.dh = 1e-8
-            self.rtol = 1e-4
+            self.rtol = 1e-3
         self.dtype = TACS.dtype
 
         self.atol = np.clip(1e-5 * self.rtol, 1e-14, 1e-8)
@@ -46,7 +46,7 @@ def setUp(self):
         )
 
         # Create stiffness (need class)
-        self.con = constitutive.IsoShellConstitutive(self.props, t=1.0, tNum=0)
+        self.con = constitutive.IsoShellConstitutive(self.props, t=1.0, tNum=0, tOffset=-0.5)
 
         # Seed random number generator in tacs for consistent test results
         elements.SeedRandomGenerator(0)
diff --git a/tests/constitutive_tests/test_smeared_composite_shell_constitutive.py b/tests/constitutive_tests/test_smeared_composite_shell_constitutive.py
index 095680992..b5983cd3d 100644
--- a/tests/constitutive_tests/test_smeared_composite_shell_constitutive.py
+++ b/tests/constitutive_tests/test_smeared_composite_shell_constitutive.py
@@ -92,6 +92,8 @@ def setUp(self):
         self.ply_fractions = np.array([0.333, 0.333, 0.333], dtype=self.dtype)
         self.ply_fraction_dv_nums = np.array([1, 2, 3], dtype=np.intc)
         self.dvs = np.array([0.1, 0.333, 0.333, 0.333], dtype=self.dtype)
+        # Distance between mid-plane and reference plane (nodes) of shell
+        self.t_offset = 0.5  # This places the nodes at the top ply
 
         # Seed random number generator in tacs for consistent test results
         elements.SeedRandomGenerator(0)
@@ -107,6 +109,7 @@ def test_constitutive_density(self):
                     self.ply_fractions,
                     self.thickness_dv_num,
                     self.ply_fraction_dv_nums,
+                    t_offset=self.t_offset,
                 )
                 fail = constitutive.TestConstitutiveDensity(
                     con,
@@ -132,6 +135,7 @@ def test_constitutive_specific_heat(self):
                     self.ply_fractions,
                     self.thickness_dv_num,
                     self.ply_fraction_dv_nums,
+                    t_offset=self.t_offset,
                 )
                 fail = constitutive.TestConstitutiveSpecificHeat(
                     con,
@@ -157,6 +161,7 @@ def test_constitutive_heat_flux(self):
                     self.ply_fractions,
                     self.thickness_dv_num,
                     self.ply_fraction_dv_nums,
+                    t_offset=self.t_offset,
                 )
                 fail = constitutive.TestConstitutiveHeatFlux(
                     con,
@@ -182,6 +187,7 @@ def test_constitutive_stress(self):
                     self.ply_fractions,
                     self.thickness_dv_num,
                     self.ply_fraction_dv_nums,
+                    t_offset=self.t_offset,
                 )
                 fail = constitutive.TestConstitutiveStress(
                     con,
@@ -207,6 +213,7 @@ def test_constitutive_thermal_strain(self):
                     self.ply_fractions,
                     self.thickness_dv_num,
                     self.ply_fraction_dv_nums,
+                    t_offset=self.t_offset,
                 )
                 fail = constitutive.TestConstitutiveThermalStrain(
                     con,
@@ -232,6 +239,7 @@ def test_constitutive_failure(self):
                     self.ply_fractions,
                     self.thickness_dv_num,
                     self.ply_fraction_dv_nums,
+                    t_offset=self.t_offset,
                 )
                 fail = constitutive.TestConstitutiveFailure(
                     con,
@@ -256,6 +264,7 @@ def test_constitutive_failure_strain_sens(self):
                     self.ply_fractions,
                     self.thickness_dv_num,
                     self.ply_fraction_dv_nums,
+                    t_offset=self.t_offset,
                 )
                 fail = constitutive.TestConstitutiveFailureStrainSens(
                     con,
diff --git a/tests/integration_tests/test_beam_bend_coupling.py b/tests/integration_tests/test_beam_bend_coupling.py
index 807134b6f..a1623fa2e 100644
--- a/tests/integration_tests/test_beam_bend_coupling.py
+++ b/tests/integration_tests/test_beam_bend_coupling.py
@@ -10,7 +10,7 @@
     Iz = 0.2
     Iy = 0.3
     J = 0.4
-    Iyz = 0.1
+    Iyz = -0.1
 Because Iyz =/= 0.0, we expect some coupling to show up in y and z bending. 
 We apply apply various tip loads test KSDisplacement, StructuralMass, MomentOfInertia, 
 and Compliance functions and sensitivities.
@@ -30,7 +30,7 @@ class ProblemTest(PyTACSTestCase.PyTACSTest):
         "z-shear_I_xy": 0.0,
         "z-shear_I_xz": 0.0,
         "z-shear_I_yy": 0.8325,
-        "z-shear_I_yz": 0.27,
+        "z-shear_I_yz": -0.27,
         "z-shear_I_zz": 0.5625,
         "z-shear_compliance": 1.346,
         "z-shear_mass": 0.27,
@@ -41,7 +41,7 @@ class ProblemTest(PyTACSTestCase.PyTACSTest):
         "y-shear_I_xy": 0.0,
         "y-shear_I_xz": 0.0,
         "y-shear_I_yy": 0.8325,
-        "y-shear_I_yz": 0.27,
+        "y-shear_I_yz": -0.27,
         "y-shear_I_zz": 0.5625,
         "y-shear_compliance": 2.006,
         "y-shear_mass": 0.27,
@@ -52,7 +52,7 @@ class ProblemTest(PyTACSTestCase.PyTACSTest):
         "x-axial_I_xy": 0.0,
         "x-axial_I_xz": 0.0,
         "x-axial_I_yy": 0.8325,
-        "x-axial_I_yz": 0.27,
+        "x-axial_I_yz": -0.27,
         "x-axial_I_zz": 0.5625,
         "x-axial_compliance": 10.0,
         "x-axial_mass": 0.27,
@@ -63,7 +63,7 @@ class ProblemTest(PyTACSTestCase.PyTACSTest):
         "x-torsion_I_xy": 0.0,
         "x-torsion_I_xz": 0.0,
         "x-torsion_I_yy": 0.8325,
-        "x-torsion_I_yz": 0.27,
+        "x-torsion_I_yz": -0.27,
         "x-torsion_I_zz": 0.5625,
         "x-torsion_compliance": 6.5,
         "x-torsion_mass": 0.27,
diff --git a/tests/integration_tests/test_rectangle_beam_offset.py b/tests/integration_tests/test_rectangle_beam_offset.py
new file mode 100644
index 000000000..ce3a63dbf
--- /dev/null
+++ b/tests/integration_tests/test_rectangle_beam_offset.py
@@ -0,0 +1,206 @@
+import os
+
+import numpy as np
+
+from pytacs_analysis_base_test import PyTACSTestCase
+from tacs import pytacs, elements, constitutive, functions
+
+"""
+This is the same test cases as `test_rectangle_beam_tractions.py`, but the beam element is offset 
+in the y and z directions, so that the beam axis no longer aligns with the nodes of the model. 
+This test ensures that the beam solution is invariant under trivial transformation: 
+
+6 noded beam model 1 meter long in x' direction.
+The cross-section is a solid rectangle with the following properties:
+    w = 0.1
+    t = 0.05
+We apply a distributed gravity load and a centrifugal load.
+We test KSDisplacement, KSFailure, StructuralMass, CenterOfMass, MomentOfInertia, and Compliance 
+functions and sensitivities.
+"""
+
+base_dir = os.path.dirname(os.path.abspath(__file__))
+bdf_file = os.path.join(base_dir, "./input_files/beam_model.bdf")
+
+ksweight = 10.0
+
+
+class ProblemTest(PyTACSTestCase.PyTACSTest):
+    N_PROCS = 2  # this is how many MPI processes to use for this TestCase.
+
+    FUNC_REFS = {
+        "gravity_I_xx": 0.0,
+        "gravity_I_xy": 0.0,
+        "gravity_I_xz": 0.0,
+        "gravity_I_yy": 1.13625,
+        "gravity_I_yz": 0.0,
+        "gravity_I_zz": 1.1278125,
+        "gravity_cgx": 0.5,
+        "gravity_cgy": -0.025,
+        "gravity_cgz": 0.05,
+        "gravity_compliance": 2000.6857974986858,
+        "gravity_ks_vmfailure": 219.37317664566595,
+        "gravity_mass": 13.5,
+        "gravity_x_disp": -0.7780125287429587,
+        "gravity_y_disp": 656.1378755069089,
+        "gravity_z_disp": 275.45079293282816,
+        "centrifugal_I_xx": 0.0,
+        "centrifugal_I_xy": 0.0,
+        "centrifugal_I_xz": 0.0,
+        "centrifugal_I_yy": 1.13625,
+        "centrifugal_I_yz": 0.0,
+        "centrifugal_I_zz": 1.1278125,
+        "centrifugal_cgx": 0.5,
+        "centrifugal_cgy": -0.025,
+        "centrifugal_cgz": 0.05,
+        "centrifugal_compliance": 9718.610964493391,
+        "centrifugal_ks_vmfailure": 402.6071884894102,
+        "centrifugal_mass": 13.5,
+        "centrifugal_x_disp": 10.454782137288209,
+        "centrifugal_y_disp": -18.840495836112435,
+        "centrifugal_z_disp": -8.065586771046936,
+    }
+
+    def setup_tacs_problems(self, comm):
+        """
+        Setup pytacs object for problems we will be testing.
+        """
+
+        # Overwrite default check values
+        if self.dtype == complex:
+            self.rtol = 1e-6
+            self.atol = 1e-6
+            self.dh = 1e-50
+        else:
+            self.rtol = 2e-1
+            self.atol = 1e-3
+            self.dh = 1e-6
+
+        # Material properties
+        rho = 2700.0  # density kg/m^3
+        E = 70.0e3  # Young's modulus (Pa)
+        nu = 0.3  # Poisson's ratio
+        ys = 2.7e3  # yield stress
+
+        # Beam dimensions
+        t = 0.05  # m
+        w = 0.1  # m
+
+        # Callback function used to setup TACS element objects and DVs
+        def elem_call_back(
+            dv_num, comp_id, comp_descript, elem_descripts, global_dvs, **kwargs
+        ):
+            # Setup (isotropic) property and constitutive objects
+            prop = constitutive.MaterialProperties(rho=rho, E=E, nu=nu, ys=ys)
+            con = constitutive.IsoRectangleBeamConstitutive(
+                prop,
+                t=t,
+                tNum=dv_num,
+                w=w,
+                wNum=dv_num + 1,
+                tOffset=0.5,
+                wOffset=-0.5,
+            )
+            refAxis = np.array([0.0, 1.0, 0.0])
+            transform = elements.BeamRefAxisTransform(refAxis)
+            # pass back the appropriate tacs element object
+            elem = elements.Beam2(transform, con)
+            return elem
+
+        # Instantiate FEA Assembler
+        struct_options = {}
+
+        fea_assembler = pytacs.pyTACS(bdf_file, comm, options=struct_options)
+
+        # Set up constitutive objects and elements
+        fea_assembler.initialize(elem_call_back)
+
+        grav_prob = fea_assembler.createStaticProblem("gravity")
+        grav_prob.addInertialLoad([-10.0, 3.0, 5.0])
+
+        rot_prob = fea_assembler.createStaticProblem("centrifugal")
+        rot_prob.addCentrifugalLoad([10.0, 3.0, 5.0], [0.5, -0.025, 0.05])
+
+        probs = [grav_prob, rot_prob]
+
+        # Set convergence to be tight for test
+        for problem in probs:
+            problem.setOption("L2Convergence", 1e-15)
+            problem.setOption("L2ConvergenceRel", 1e-15)
+
+        # Add Functions
+        for problem in probs:
+            problem.addFunction("mass", functions.StructuralMass)
+            problem.addFunction("compliance", functions.Compliance)
+            problem.addFunction("ks_vmfailure", functions.KSFailure, ksWeight=ksweight)
+            problem.addFunction(
+                "x_disp",
+                functions.KSDisplacement,
+                ksWeight=ksweight,
+                direction=[10.0, 0.0, 0.0],
+            )
+            problem.addFunction(
+                "y_disp",
+                functions.KSDisplacement,
+                ksWeight=ksweight,
+                direction=[0.0, 10.0, 0.0],
+            )
+            problem.addFunction(
+                "z_disp",
+                functions.KSDisplacement,
+                ksWeight=ksweight,
+                direction=[0.0, 0.0, 10.0],
+            )
+            problem.addFunction(
+                "cgx", functions.CenterOfMass, direction=[1.0, 0.0, 0.0]
+            )
+            problem.addFunction(
+                "cgy", functions.CenterOfMass, direction=[0.0, 1.0, 0.0]
+            )
+            problem.addFunction(
+                "cgz", functions.CenterOfMass, direction=[0.0, 0.0, 1.0]
+            )
+            problem.addFunction(
+                "I_xx",
+                functions.MomentOfInertia,
+                direction1=[1.0, 0.0, 0.0],
+                direction2=[1.0, 0.0, 0.0],
+                aboutCM=True,
+            )
+            problem.addFunction(
+                "I_xy",
+                functions.MomentOfInertia,
+                direction1=[1.0, 0.0, 0.0],
+                direction2=[0.0, 1.0, 0.0],
+                aboutCM=True,
+            )
+            problem.addFunction(
+                "I_xz",
+                functions.MomentOfInertia,
+                direction1=[1.0, 0.0, 0.0],
+                direction2=[0.0, 0.0, 1.0],
+                aboutCM=True,
+            )
+            problem.addFunction(
+                "I_yy",
+                functions.MomentOfInertia,
+                direction1=[0.0, 1.0, 0.0],
+                direction2=[0.0, 1.0, 0.0],
+                aboutCM=True,
+            )
+            problem.addFunction(
+                "I_yz",
+                functions.MomentOfInertia,
+                direction1=[0.0, 1.0, 0.0],
+                direction2=[0.0, 0.0, 1.0],
+                aboutCM=True,
+            )
+            problem.addFunction(
+                "I_zz",
+                functions.MomentOfInertia,
+                direction1=[0.0, 0.0, 1.0],
+                direction2=[0.0, 0.0, 1.0],
+                aboutCM=True,
+            )
+
+        return probs, fea_assembler
diff --git a/tests/integration_tests/test_shell_blade_stiffened_plate_quad.py b/tests/integration_tests/test_shell_blade_stiffened_plate_quad.py
index fda8da484..463d7d3bb 100644
--- a/tests/integration_tests/test_shell_blade_stiffened_plate_quad.py
+++ b/tests/integration_tests/test_shell_blade_stiffened_plate_quad.py
@@ -32,39 +32,39 @@ class ProblemTest(PyTACSTestCase.PyTACSTest):
     N_PROCS = 2  # this is how many MPI processes to use for this TestCase.
 
     FUNC_REFS = {
-        "point_load_Ixx": 1.4591536458333065,
-        "point_load_Ixy": 3.907985046680551e-14,
+        "point_load_Ixx": 1.4589304315475928,
+        "point_load_Ixy": 3.8191672047105385e-14,
         "point_load_Ixz": 0.0,
-        "point_load_Iyy": 1.459153645833294,
+        "point_load_Iyy": 1.4589304315475786,
         "point_load_Iyz": 0.0,
         "point_load_Izz": 2.916666666666588,
         "point_load_cgx": 0.500000000000004,
         "point_load_cgy": 0.500000000000004,
-        "point_load_cgz": 0.0,
+        "point_load_cgz": -0.0035714285714285718,
         "point_load_compliance": 71.87161795633577,
         "point_load_ks_vmfailure": 0.3304517634314118,
         "point_load_mass": 17.5,
-        "pressure_Ixx": 1.4591536458333065,
+        "pressure_Ixx": 1.4589304315475928,
         "pressure_Ixy": 3.907985046680551e-14,
         "pressure_Ixz": 0.0,
-        "pressure_Iyy": 1.459153645833294,
+        "pressure_Iyy": 1.4589304315475786,
         "pressure_Iyz": 0.0,
         "pressure_Izz": 2.916666666666588,
         "pressure_cgx": 0.500000000000004,
         "pressure_cgy": 0.500000000000004,
-        "pressure_cgz": 0.0,
+        "pressure_cgz": -0.0035714285714285718,
         "pressure_compliance": 377.11604579572935,
         "pressure_ks_vmfailure": 0.9085085771277308,
         "pressure_mass": 17.5,
-        "gravity_Ixx": 1.4591536458333065,
+        "gravity_Ixx": 1.4589304315475928,
         "gravity_Ixy": 3.907985046680551e-14,
         "gravity_Ixz": 0.0,
-        "gravity_Iyy": 1.459153645833294,
+        "gravity_Iyy": 1.4589304315475786,
         "gravity_Iyz": 0.0,
         "gravity_Izz": 2.916666666666588,
         "gravity_cgx": 0.500000000000004,
         "gravity_cgy": 0.500000000000004,
-        "gravity_cgz": 0.0,
+        "gravity_cgz": -0.0035714285714285718,
         "gravity_compliance": 11.114479357783475,
         "gravity_ks_vmfailure": 0.0908244150089233,
         "gravity_mass": 17.5,
diff --git a/tests/integration_tests/test_shell_comp_offset.py b/tests/integration_tests/test_shell_comp_offset.py
new file mode 100644
index 000000000..1f50bc48b
--- /dev/null
+++ b/tests/integration_tests/test_shell_comp_offset.py
@@ -0,0 +1,254 @@
+import os
+
+import numpy as np
+
+from pytacs_analysis_base_test import PyTACSTestCase
+from tacs import pytacs, elements, constitutive, functions
+
+"""
+This is the same geometry as `test_shell_plate_quad.py`, but the plate is offset 
+in the z direction, so that the shell plane no longer aligns with the nodes of the model.
+Tests a smeared laminate shell model with the following layup: [0, 45, 30].
+Agravity load in the x direction is applied.
+tests KSDisplacement, KSFailure, StructuralMass, CenterOfMass, MomentOfInertia, and Compliance functions
+and sensitivities.
+We also test a ply fraction summation constraint using the DVConstraint class.
+"""
+
+base_dir = os.path.dirname(os.path.abspath(__file__))
+bdf_file = os.path.join(base_dir, "./input_files/comp_plate.bdf")
+
+ksweight = 10.0
+
+
+class ProblemTest(PyTACSTestCase.PyTACSTest):
+    N_PROCS = 2  # this is how many MPI processes to use for this TestCase.
+
+    FUNC_REFS = {
+        "gravity_I_xx": 0.38750005449218716,
+        "gravity_I_xy": 5.551115123125783e-17,
+        "gravity_I_xz": -2.710505431213761e-20,
+        "gravity_I_yy": 0.02421880449218755,
+        "gravity_I_yz": -2.168404344971009e-19,
+        "gravity_I_zz": 0.4117187499999999,
+        "gravity_cgx": 0.25000000000000006,
+        "gravity_cgy": 1.0000000000000002,
+        "gravity_cgz": -0.00037500000000000006,
+        "gravity_compliance": 3.4205235434811567e-06,
+        "gravity_ks_TsaiWufailure": 0.2259320352899525,
+        "gravity_mass": 1.1624999999999999,
+        "gravity_x_disp": 1.625976981601372e-06,
+        "gravity_y_disp": -1.2736463024905141e-07,
+        "gravity_z_disp": -0.000920864663486706,
+        "centrifugal_I_xx": 0.38750005449218716,
+        "centrifugal_I_xy": 5.551115123125783e-17,
+        "centrifugal_I_xz": -2.710505431213761e-20,
+        "centrifugal_I_yy": 0.02421880449218755,
+        "centrifugal_I_yz": -2.168404344971009e-19,
+        "centrifugal_I_zz": 0.4117187499999999,
+        "centrifugal_cgx": 0.25000000000000006,
+        "centrifugal_cgy": 1.0000000000000002,
+        "centrifugal_cgz": -0.00037500000000000006,
+        "centrifugal_compliance": 0.23430251098271168,
+        "centrifugal_ks_TsaiWufailure": 0.24632761512178053,
+        "centrifugal_mass": 1.1624999999999999,
+        "centrifugal_x_disp": 0.0001570586015183079,
+        "centrifugal_y_disp": 8.748943337216315e-05,
+        "centrifugal_z_disp": 0.0966983775511979,
+    }
+
+    def setup_tacs_problems(self, comm):
+        """
+        Setup pytacs object for problems we will be testing.
+        """
+
+        # Overwrite default check values
+        if self.dtype == complex:
+            self.rtol = 1e-6
+            self.atol = 1e-6
+            self.dh = 1e-50
+        else:
+            self.rtol = 1e-3
+            self.atol = 1e-3
+            self.dh = 1e-7
+
+        fea_assembler = pytacs.pyTACS(bdf_file, comm)
+
+        # Material properties
+        rho = 1550.0
+        specific_heat = 921.096
+        E1 = 54e9
+        E2 = 18e9
+        nu12 = 0.25
+        G12 = 9e9
+        G13 = 9e9
+        Xt = 2410.0e6
+        Xc = 1040.0e6
+        Yt = 73.0e6
+        Yc = 173.0e6
+        S12 = 71.0e6
+        cte = 24.0e-6
+        kappa = 230.0
+        ply_thickness = 1.25e-3
+
+        ortho_prop = constitutive.MaterialProperties(
+            rho=rho,
+            specific_heat=specific_heat,
+            E1=E1,
+            E2=E2,
+            nu12=nu12,
+            G12=G12,
+            G13=G13,
+            G23=G13,
+            Xt=Xt,
+            Xc=Xc,
+            Yt=Yt,
+            Yc=Yc,
+            S12=S12,
+            alpha=cte,
+            kappa=kappa,
+        )
+        ortho_ply = constitutive.OrthotropicPly(ply_thickness, ortho_prop)
+
+        # Shell thickness
+        t = 7.5e-4  # m
+
+        # Define reference axis for local shell stresses
+        refAxis = np.array([1.0, 0.0, 0.0])
+        transform = elements.ShellRefAxisTransform(refAxis)
+
+        # Callback function used to setup TACS element objects and DVs
+        def elemCallBack(
+            dvNum, compID, compDescript, elemDescripts, specialDVs, **kwargs
+        ):
+            # Create the orthotropic layup
+            ortho_layup = [ortho_ply] * 3
+            ply_angles = np.deg2rad([0.0, 45.0, 30.0]).astype(self.dtype)
+            ply_fractions = np.array([1.0 / 3.0] * 3, dtype=self.dtype)
+            thickness_dv_num = dvNum
+            ply_fraction_dv_nums = np.array(
+                [dvNum + 1, dvNum + 2, dvNum + 3], dtype=np.intc
+            )
+
+            con = constitutive.SmearedCompositeShellConstitutive(
+                ortho_layup,
+                t,
+                ply_angles,
+                ply_fractions,
+                thickness_dv_num=thickness_dv_num,
+                ply_fraction_dv_nums=ply_fraction_dv_nums,
+                t_offset=0.5,
+            )
+
+            # For each element type in this component,
+            # pass back the appropriate tacs element object
+            elemList = []
+            for elemDescript in elemDescripts:
+                if elemDescript in ["CQUAD4", "CQUADR"]:
+                    elem = elements.Quad4Shell(transform, con)
+                elif elemDescript in ["CTRIA3", "CTRIAR"]:
+                    elem = elements.Tri3Shell(transform, con)
+                else:
+                    print("Uh oh, '%s' not recognized" % (elemDescript))
+                elemList.append(elem)
+
+            return elemList
+
+        # Set up constitutive objects and elements
+        fea_assembler.initialize(elemCallBack)
+
+        # Read in forces from BDF and create tacs struct problems
+        grav_prob = fea_assembler.createStaticProblem("gravity")
+        grav_prob.addInertialLoad(1000 * [9.81, 0.0, 0.0])
+
+        rot_prob = fea_assembler.createStaticProblem("centrifugal")
+        rot_prob.addCentrifugalLoad([10.0, 3.0, 5.0], [0.25, 1.0, -0.000375])
+
+        probs = [grav_prob, rot_prob]
+
+        # Set convergence to be tight for test
+        for problem in probs:
+            problem.setOption("L2Convergence", 1e-15)
+            problem.setOption("L2ConvergenceRel", 1e-15)
+
+        # Add Functions
+        for problem in probs:
+            problem.addFunction("mass", functions.StructuralMass)
+            problem.addFunction("compliance", functions.Compliance)
+            problem.addFunction(
+                "ks_TsaiWufailure",
+                functions.KSFailure,
+                ksWeight=ksweight,
+                ftype="discrete",
+            )
+            problem.addFunction(
+                "x_disp",
+                functions.KSDisplacement,
+                ksWeight=ksweight,
+                direction=[10.0, 0.0, 0.0],
+            )
+            problem.addFunction(
+                "y_disp",
+                functions.KSDisplacement,
+                ksWeight=ksweight,
+                direction=[0.0, 10.0, 0.0],
+            )
+            problem.addFunction(
+                "z_disp",
+                functions.KSDisplacement,
+                ksWeight=ksweight,
+                direction=[0.0, 0.0, 10.0],
+            )
+            problem.addFunction(
+                "cgx", functions.CenterOfMass, direction=[1.0, 0.0, 0.0]
+            )
+            problem.addFunction(
+                "cgy", functions.CenterOfMass, direction=[0.0, 1.0, 0.0]
+            )
+            problem.addFunction(
+                "cgz", functions.CenterOfMass, direction=[0.0, 0.0, 1.0]
+            )
+            problem.addFunction(
+                "I_xx",
+                functions.MomentOfInertia,
+                direction1=[1.0, 0.0, 0.0],
+                direction2=[1.0, 0.0, 0.0],
+                aboutCM=True,
+            )
+            problem.addFunction(
+                "I_xy",
+                functions.MomentOfInertia,
+                direction1=[1.0, 0.0, 0.0],
+                direction2=[0.0, 1.0, 0.0],
+                aboutCM=True,
+            )
+            problem.addFunction(
+                "I_xz",
+                functions.MomentOfInertia,
+                direction1=[1.0, 0.0, 0.0],
+                direction2=[0.0, 0.0, 1.0],
+                aboutCM=True,
+            )
+            problem.addFunction(
+                "I_yy",
+                functions.MomentOfInertia,
+                direction1=[0.0, 1.0, 0.0],
+                direction2=[0.0, 1.0, 0.0],
+                aboutCM=True,
+            )
+            problem.addFunction(
+                "I_yz",
+                functions.MomentOfInertia,
+                direction1=[0.0, 1.0, 0.0],
+                direction2=[0.0, 0.0, 1.0],
+                aboutCM=True,
+            )
+            problem.addFunction(
+                "I_zz",
+                functions.MomentOfInertia,
+                direction1=[0.0, 0.0, 1.0],
+                direction2=[0.0, 0.0, 1.0],
+                aboutCM=True,
+            )
+
+        return probs, fea_assembler
diff --git a/tests/integration_tests/test_shell_comp_smeared.py b/tests/integration_tests/test_shell_comp_smeared.py
index ba187f0e9..b6b2a70aa 100644
--- a/tests/integration_tests/test_shell_comp_smeared.py
+++ b/tests/integration_tests/test_shell_comp_smeared.py
@@ -10,7 +10,7 @@
 Two load cases are tested: an in-plane tension and out-of-plane shear.
 This test is identical to test_shell_comp_unbalanced.py except since the laminate 
 is smeared all stacking sequence dependence is neglected.
-tests KSDisplacement, KSFailure, StructuralMass, and Compliance functions
+tests KSDisplacement, KSFailure, StructuralMass, CenterOfMass, MomentOfInertia, and Compliance functions
 and sensitivities.
 We also test a ply fraction summation constraint using the DVConstraint class.
 """
@@ -25,19 +25,36 @@ class ProblemTest(PyTACSTestCase.PyTACSTest):
     N_PROCS = 2  # this is how many MPI processes to use for this TestCase.
 
     FUNC_REFS = {
-        "Tension_compliance": 8433.692626389691,
-        "Tension_ks_TsaiWufailure": 4.352340467252131,
+        "Tension_I_xx": 0.3875000544921874,
+        "Tension_I_xy": 5.551115123125783e-17,
+        "Tension_I_xz": 0.0,
+        "Tension_I_yy": 0.02421880449218751,
+        "Tension_I_yz": 0.0,
+        "Tension_I_zz": 0.4117187499999999,
+        "Tension_cgx": 0.25000000000000006,
+        "Tension_cgy": 1.0000000000000002,
+        "Tension_cgz": 0.0,
+        "Tension_compliance": 8433.69262638969,
+        "Tension_ks_TsaiWufailure": 4.35234046725213,
         "Tension_mass": 1.1624999999999999,
         "Tension_x_disp": 0.04604283873354243,
-        "Tension_y_disp": -0.01468482384786706,
+        "Tension_y_disp": -0.014684823847867035,
         "Tension_z_disp": 0.0,
-        "VertShear_compliance": 0.0001081728488799184,
-        "VertShear_ks_TsaiWufailure": 0.22626387488407798,
+        "VertShear_I_xx": 0.3875000544921874,
+        "VertShear_I_xy": 5.551115123125783e-17,
+        "VertShear_I_xz": 0.0,
+        "VertShear_I_yy": 0.02421880449218751,
+        "VertShear_I_yz": 0.0,
+        "VertShear_I_zz": 0.4117187499999999,
+        "VertShear_cgx": 0.25000000000000006,
+        "VertShear_cgy": 1.0000000000000002,
+        "VertShear_cgz": 0.0,
+        "VertShear_compliance": 0.00010817284888043632,
+        "VertShear_ks_TsaiWufailure": 0.22626387488408603,
         "VertShear_mass": 1.1624999999999999,
         "VertShear_x_disp": 0.0,
         "VertShear_y_disp": 0.0,
-        "VertShear_z_disp": 0.0054598659927654735,
-        "ply_fractions_sum": 1.0,
+        "VertShear_z_disp": 0.0054598659927803965,
     }
 
     def setup_tacs_problems(self, comm):
@@ -176,15 +193,58 @@ def elemCallBack(
                 ksWeight=ksweight,
                 direction=[0.0, 0.0, 10.0],
             )
+            problem.addFunction(
+                "cgx", functions.CenterOfMass, direction=[1.0, 0.0, 0.0]
+            )
+            problem.addFunction(
+                "cgy", functions.CenterOfMass, direction=[0.0, 1.0, 0.0]
+            )
+            problem.addFunction(
+                "cgz", functions.CenterOfMass, direction=[0.0, 0.0, 1.0]
+            )
+            problem.addFunction(
+                "I_xx",
+                functions.MomentOfInertia,
+                direction1=[1.0, 0.0, 0.0],
+                direction2=[1.0, 0.0, 0.0],
+                aboutCM=True,
+            )
+            problem.addFunction(
+                "I_xy",
+                functions.MomentOfInertia,
+                direction1=[1.0, 0.0, 0.0],
+                direction2=[0.0, 1.0, 0.0],
+                aboutCM=True,
+            )
+            problem.addFunction(
+                "I_xz",
+                functions.MomentOfInertia,
+                direction1=[1.0, 0.0, 0.0],
+                direction2=[0.0, 0.0, 1.0],
+                aboutCM=True,
+            )
+            problem.addFunction(
+                "I_yy",
+                functions.MomentOfInertia,
+                direction1=[0.0, 1.0, 0.0],
+                direction2=[0.0, 1.0, 0.0],
+                aboutCM=True,
+            )
+            problem.addFunction(
+                "I_yz",
+                functions.MomentOfInertia,
+                direction1=[0.0, 1.0, 0.0],
+                direction2=[0.0, 0.0, 1.0],
+                aboutCM=True,
+            )
+            problem.addFunction(
+                "I_zz",
+                functions.MomentOfInertia,
+                direction1=[0.0, 0.0, 1.0],
+                direction2=[0.0, 0.0, 1.0],
+                aboutCM=True,
+            )
 
         tacs_probs = list(tacs_probs)
 
-        # Add linear constraint on ply fraction summation
-        constr = fea_assembler.createDVConstraint("ply_fractions")
-        allComponents = fea_assembler.selectCompIDs()
-        constr.addConstraint(
-            "sum", allComponents, dvIndices=[1, 2, 3], dvWeights=[1.0, 1.0, 1.0]
-        )
-        tacs_probs.append(constr)
-
         return tacs_probs, fea_assembler
diff --git a/tests/integration_tests/test_shell_plate_offset.py b/tests/integration_tests/test_shell_plate_offset.py
new file mode 100644
index 000000000..7e1bb0676
--- /dev/null
+++ b/tests/integration_tests/test_shell_plate_offset.py
@@ -0,0 +1,154 @@
+import os
+
+import numpy as np
+
+from pytacs_analysis_base_test import PyTACSTestCase
+from tacs import pytacs, elements, constitutive, functions
+
+"""
+This is the same test cases as `test_shell_plate_quad.py`, but the plate is offset 
+in the z direction, so that the shell plane no longer aligns with the nodes of the model.
+The nominal case is a 1m x 1m flat plate under a 100G gravity load. The
+perimeter of the plate is fixed in all 6 degrees of freedom. The plate comprises
+100 CQUAD4 elements and test KSFailure, StructuralMass, CenterOfMass, MomentOfInertia,
+and Compliance functions and sensitivities.
+"""
+
+base_dir = os.path.dirname(os.path.abspath(__file__))
+bdf_file = os.path.join(base_dir, "./input_files/plate.bdf")
+
+# KS function weight
+ksweight = 10.0
+
+
+class ProblemTest(PyTACSTestCase.PyTACSTest):
+    N_PROCS = 2  # this is how many MPI processes to use for this TestCase.
+
+    FUNC_REFS = {
+        "gravity_Ixx": 1.041692708333326,
+        "gravity_Ixy": 4.884981308350689e-15,
+        "gravity_Ixz": 0.0,
+        "gravity_Iyy": 1.0416927083333283,
+        "gravity_Iyz": 0.0,
+        "gravity_Izz": 2.0833333333333233,
+        "gravity_cgx": 0.5,
+        "gravity_cgy": 0.5,
+        "gravity_cgz": 0.0025,
+        "gravity_compliance": 67.59168942441727,
+        "gravity_ks_vmfailure": 0.11479343495254533,
+        "gravity_mass": 12.5,
+    }
+
+    def setup_tacs_problems(self, comm):
+        """
+        Setup pytacs object for problems we will be testing.
+        """
+
+        # Overwrite default check values
+        if self.dtype == complex:
+            self.rtol = 1e-8
+            self.atol = 1e-8
+            self.dh = 1e-50
+        else:
+            self.rtol = 2e-1
+            self.atol = 1e-4
+            self.dh = 1e-6
+
+        # Instantiate FEA Assembler
+        struct_options = {}
+
+        fea_assembler = pytacs.pyTACS(bdf_file, comm, options=struct_options)
+
+        def elem_call_back(
+            dv_num, comp_id, comp_descript, elem_descripts, global_dvs, **kwargs
+        ):
+            # Material properties
+            rho = 2500.0  # density kg/m^3
+            E = 70e9  # Young's modulus (Pa)
+            nu = 0.3  # Poisson's ratio
+            ys = 464.0e6  # yield stress
+
+            # Plate geometry
+            tplate = 0.005  # 5 mm
+
+            # Set up property model
+            prop = constitutive.MaterialProperties(rho=rho, E=E, nu=nu, ys=ys)
+            # Set up constitutive model
+            con = constitutive.IsoShellConstitutive(
+                prop, t=tplate, tNum=dv_num, tOffset=-0.5
+            )
+            transform = None
+            # Set up element
+            elem = elements.Quad4Shell(transform, con)
+            scale = [100.0]
+            return elem, scale
+
+        # Set up constitutive objects and elements
+        fea_assembler.initialize(elem_call_back)
+
+        tacs_probs = []
+
+        # Add pressure to entire plate
+        sp = fea_assembler.createStaticProblem(name="gravity")
+        g = np.array([0.0, 0.0, -981.0], dtype=self.dtype)
+        sp.addInertialLoad(g)
+        tacs_probs.append(sp)
+
+        # Add Functions
+        for problem in tacs_probs:
+            problem.addFunction("mass", functions.StructuralMass)
+            problem.addFunction("ks_vmfailure", functions.KSFailure, ksWeight=ksweight)
+            problem.addFunction("compliance", functions.Compliance)
+            problem.addFunction(
+                "cgx", functions.CenterOfMass, direction=[1.0, 0.0, 0.0]
+            )
+            problem.addFunction(
+                "cgy", functions.CenterOfMass, direction=[0.0, 1.0, 0.0]
+            )
+            problem.addFunction(
+                "cgz", functions.CenterOfMass, direction=[0.0, 0.0, 1.0]
+            )
+            problem.addFunction(
+                "Ixx",
+                functions.MomentOfInertia,
+                direction1=[1.0, 0.0, 0.0],
+                direction2=[1.0, 0.0, 0.0],
+                aboutCM=True,
+            )
+            problem.addFunction(
+                "Ixy",
+                functions.MomentOfInertia,
+                direction1=[1.0, 0.0, 0.0],
+                direction2=[0.0, 1.0, 0.0],
+                aboutCM=True,
+            )
+            problem.addFunction(
+                "Ixz",
+                functions.MomentOfInertia,
+                direction1=[1.0, 0.0, 0.0],
+                direction2=[0.0, 0.0, 1.0],
+                aboutCM=True,
+            )
+            problem.addFunction(
+                "Iyy",
+                functions.MomentOfInertia,
+                direction1=[0.0, 1.0, 0.0],
+                direction2=[0.0, 1.0, 0.0],
+                aboutCM=True,
+            )
+            problem.addFunction(
+                "Iyz",
+                functions.MomentOfInertia,
+                direction1=[0.0, 0.0, 1.0],
+                direction2=[0.0, 1.0, 0.0],
+                aboutCM=True,
+            )
+            problem.addFunction(
+                "Izz",
+                functions.MomentOfInertia,
+                direction1=[0.0, 0.0, 1.0],
+                direction2=[0.0, 0.0, 1.0],
+                aboutCM=True,
+            )
+
+        return tacs_probs, fea_assembler
diff --git a/tests/integration_tests/test_trans_beam.py b/tests/integration_tests/test_trans_beam.py
index 3e436fcb6..b638022f4 100644
--- a/tests/integration_tests/test_trans_beam.py
+++ b/tests/integration_tests/test_trans_beam.py
@@ -23,13 +23,13 @@ class ProblemTest(PyTACSTestCase.PyTACSTest):
     N_PROCS = 2  # this is how many MPI processes to use for this TestCase.
 
     FUNC_REFS = {
-        "ramp_compliance": 10.507315462331803,
+        "ramp_compliance": 10.917589919679486,
         "ramp_x_disp": 0.06931471805599457,
-        "ramp_y_disp": 12.186052999423167,
+        "ramp_y_disp": 12.191784286714391,
         "ramp_z_disp": 0.06931471805599457,
-        "sinusoid_compliance": 22.093569888841497,
+        "sinusoid_compliance": 23.113858857368683,
         "sinusoid_x_disp": 0.06931471805599457,
-        "sinusoid_y_disp": 25.654265895487846,
+        "sinusoid_y_disp": 25.47755094553053,
         "sinusoid_z_disp": 0.06931471805599457,
     }
 
@@ -49,9 +49,7 @@ def setup_tacs_problems(self, comm):
             self.dh = 1e-6
 
         # Instantiate FEA Assembler
-        struct_options = {}
-
-        fea_assembler = pytacs.pyTACS(bdf_file, comm, options=struct_options)
+        fea_assembler = pytacs.pyTACS(bdf_file, comm)
 
         # Material properties
         rho = 27.0  # density kg/m^3
@@ -78,9 +76,6 @@ def elem_call_back(
             elem = elements.Beam2(transform, con)
             return elem
 
-        # Instantiate FEA Assembler
-        fea_assembler = pytacs.pyTACS(bdf_file, comm)
-
         # Set up constitutive objects and elements
         fea_assembler.initialize(elem_call_back)
 
@@ -90,8 +85,8 @@ def elem_call_back(
         tacs_probs = tacs_probs.values()
         # Set convergence to be tight for test
         for problem in tacs_probs:
-            problem.setOption("L2Convergence", 1e-15)
-            problem.setOption("L2ConvergenceRel", 1e-15)
+            problem.setOption("L2Convergence", 1e-16)
+            problem.setOption("L2ConvergenceRel", 1e-16)
 
         # Add Functions
         for problem in tacs_probs: