diff --git a/beam.hpp b/beam.hpp
index 9557074..5380171 100644
--- a/beam.hpp
+++ b/beam.hpp
@@ -31,13 +31,13 @@ struct CylindricalBeamDimensions {
 };
 
 struct ElementMaterial {
-  float G; // poisson's ratio
-  float E; // ym in pascal
+  float poissonsRatio; // poisson's ratio
+  float E;             // ym in pascal
   ElementMaterial(const float &poissonsRatio, const float &youngsModulus)
-      : G(poissonsRatio), E(youngsModulus) {
+      : poissonsRatio(poissonsRatio), E(youngsModulus) {
     assert(poissonsRatio <= 0.5 && poissonsRatio >= -1);
   }
-  ElementMaterial() : G(0.3), E(200 * 1e9) {}
+  ElementMaterial() : poissonsRatio(0.3), E(200 * 1e9) {}
 };
 
 #endif // BEAM_HPP
diff --git a/beamformfinder.cpp b/beamformfinder.cpp
index 5aa2332..826817c 100644
--- a/beamformfinder.cpp
+++ b/beamformfinder.cpp
@@ -756,9 +756,8 @@ double FormFinder::computeTotalPotentialEnergy() const {
                            element.properties.A / (2 * element.initialLength);
     const double theta1Diff = element.rotationalDisplacements_jplus1.theta1 -
                               element.rotationalDisplacements_j.theta1;
-    const double torsionalPE = element.properties.material.G *
-                               element.properties.J * pow(theta1Diff, 2) /
-                               (2 * element.initialLength);
+    const double torsionalPE = element.properties.G * element.properties.J *
+                               pow(theta1Diff, 2) / (2 * element.initialLength);
     const double &theta2_j = element.rotationalDisplacements_j.theta2;
     const double &theta2_jplus1 = element.rotationalDisplacements_jplus1.theta2;
     const double firstBendingPE_inBracketsTerm = 4 * pow(theta2_j, 2) +
diff --git a/elementalmesh.cpp b/elementalmesh.cpp
index 3d0b839..6bed494 100644
--- a/elementalmesh.cpp
+++ b/elementalmesh.cpp
@@ -419,6 +419,7 @@ double computeAngle(const VectorType &vector0, const VectorType &vector1,
 void Element::Properties::computeMaterialProperties(
     const ElementMaterial &material) {
   this->material = material;
+  this->G = material.E / (2 * (1 + material.poissonsRatio));
 }
 
 void Element::Properties::computeDimensionsProperties(
@@ -459,7 +460,7 @@ Element::Properties::Properties(const CrossSectionType &dimensions,
 Element::Properties::Properties(const std::array<double, 6> &arr) {
   assert(arr.size() == 6);
   material.E = arr[0];
-  material.G = arr[1];
+  G = arr[1];
   A = arr[2];
   I2 = arr[3];
   I3 = arr[4];
@@ -467,14 +468,13 @@ Element::Properties::Properties(const std::array<double, 6> &arr) {
 }
 
 std::array<double, 6> Element::Properties::toArray() const {
-  return std::array<double, 6>({material.E, material.G, A,
-                                static_cast<double>(I2),
-                                static_cast<double>(I3), J});
+  return std::array<double, 6>(
+      {material.E, G, A, static_cast<double>(I2), static_cast<double>(I3), J});
 }
 
 void Element::updateRigidity() {
   rigidity.axial = properties.material.E * properties.A / initialLength;
-  rigidity.torsional = properties.material.G * properties.J / initialLength;
+  rigidity.torsional = properties.G * properties.J / initialLength;
   rigidity.firstBending =
       2 * properties.material.E * properties.I2 / initialLength;
   rigidity.secondBending =
diff --git a/elementalmesh.hpp b/elementalmesh.hpp
index 4c8a572..1436b2c 100644
--- a/elementalmesh.hpp
+++ b/elementalmesh.hpp
@@ -56,6 +56,7 @@ struct Element {
   struct Properties {
     CrossSectionType dimensions;
     ElementMaterial material;
+    double G{0};
     double A{0};  // cross sectional area
     double I2{0}; // second moment of inertia
     double I3{0}; // third moment of inertia