diff --git a/dakota_interface.py b/dakota_interface.py
index 3debb82..c9162a4 100644
--- a/dakota_interface.py
+++ b/dakota_interface.py
@@ -69,8 +69,12 @@ def process(fn_json_args):
             scrnout = interface_args['scrnout']
         ms = sga.solve(fn_main, analysis, ppcmd, solver, scrnout=scrnout)
 
+        if 'ref_center' in interface_args.keys():
+            if interface_args['ref_center'] == 'sc':
+                ms.offsetBeamRefCenter(ms.shear_center[1], ms.shear_center[2])
 
         # Calculate actual responses
+        rc = None
         for k, v in interface_args['post_process'].items():
             # print('')
             # print('v[0]:', v[0])
diff --git a/sg.py b/sg.py
index 47c1956..37a5075 100644
--- a/sg.py
+++ b/sg.py
@@ -100,7 +100,7 @@ def summary(self):
                     print(row)
 
 
-    def getBeamProperty(self, label, center=None):
+    def getBeamProperty(self, label):
         """Get beam properties using specific labels.
 
         Parameters
@@ -145,13 +145,13 @@ def getBeamProperty(self, label, center=None):
         stf_r = self.stiffness_refined
         cmp_r = self.compliance_refined
 
-        if type(center).__name__ == 'str':
-            if center == 'tc':
-                mm, stf_c, cmp_c, stf_r, cmp_r = calcOffsetBeamProperty(self.tension_center[1], self.tension_center[2])
-            elif center == 'sc':
-                mm, stf_c, cmp_c, stf_r, cmp_r = calcOffsetBeamProperty(self.shear_center[1], self.shear_center[2])
-            elif center == 'mc':
-                mm, stf_c, cmp_c, stf_r, cmp_r = calcOffsetBeamProperty(self.mass_center[1], self.mass_center[2])
+        # if type(center).__name__ == 'str':
+        #     if center == 'tc':
+        #         mm, stf_c, cmp_c, stf_r, cmp_r = calcOffsetBeamProperty(self.tension_center[1], self.tension_center[2])
+        #     elif center == 'sc':
+        #         mm, stf_c, cmp_c, stf_r, cmp_r = calcOffsetBeamProperty(self.shear_center[1], self.shear_center[2])
+        #     elif center == 'mc':
+        #         mm, stf_c, cmp_c, stf_r, cmp_r = calcOffsetBeamProperty(self.mass_center[1], self.mass_center[2])
 
         # Mass
         if label.startswith('ms'):
@@ -266,6 +266,75 @@ def calcOffsetBeamProperty(self, offset_x2=0.0, offset_x3=0.0):
         return mm_o, stf_4_trf, cmp_4_trf, stf_6_trf, cmp_6_trf
 
 
+    def offsetBeamRefCenter(self, offset_x2, offset_x3):
+        """Offset the beam reference center and recalculate beam properties.
+
+        Parameters
+        ----------
+        offset_x2 : float
+            x2 of the offset of the new center with respect to the current one.
+        offset_x3 : float
+            x3 of the offset of the new center with respect to the current one.
+
+        """
+
+        # offset = [0.0, offset[0], offset[1]]
+
+        # Offset mass matrix
+        mm_o = np.asarray(self.mass_mc)
+        if (offset_x2 != self.mass_center[1]) or (offset_x3 != self.mass_center[2]):
+            # mm_c = np.asarray(self.mass_mc)
+            mu = mm_o[0, 0]
+            mi_c = mm_o[3:, 3:]
+
+            x2 = self.mass_center[1] - offset_x2
+            x3 = self.mass_center[2] - offset_x3
+            r_tilde = np.array([
+                [0, -x3, x2],
+                [x3, 0, 0],
+                [-x2, 0, 0]
+            ])
+
+            mm_o[:3, 3:] = mu * r_tilde.T
+            mm_o[3:, :3] = mu * r_tilde
+
+            # I_o = I_c + m * r_tilde.r_tilde^T
+            mm_o[3:, 3:] = mm_o[3:, 3:] + mu * np.dot(r_tilde, r_tilde.T)
+        self.mass_origin = mm_o
+        self.mass_center[1] -= offset_x2
+        self.mass_center[2] -= offset_x3
+
+
+        # Offset stiffness and compliance
+        trfm_4 = np.eye(4)
+        trfm_6 = np.eye(6)
+
+        trfm_4[2, 0] = offset_x3
+        trfm_4[3, 0] = -offset_x2
+
+        trfm_6[4, 0] = offset_x3
+        trfm_6[5, 0] = -offset_x2
+        trfm_6[3, 1] = -offset_x3
+        trfm_6[3, 2] = offset_x2
+
+        cmp_4 = np.asarray(self.compliance)
+        cmp_6 = np.asarray(self.compliance_refined)
+
+        self.compliance = np.dot(trfm_4.T, np.dot(cmp_4, trfm_4))
+        self.compliance_refined = np.dot(trfm_6.T, np.dot(cmp_6, trfm_6))
+
+        self.stiffness = np.linalg.inv(self.compliance)
+        self.stiffness_refined = np.linalg.inv(self.compliance_refined)
+
+        self.tension_center[1] -= offset_x2
+        self.tension_center[2] -= offset_x3
+
+        self.shear_center[1] -= offset_x2
+        self.shear_center[2] -= offset_x3
+
+        return
+
+