diff --git a/examples/solidHeatTransferScript/HeatConduction2DFin/HeatConduction2DFin.html b/examples/solidHeatTransferScript/HeatConduction2DFin/HeatConduction2DFin.html
index 8c98ab1..6fb6c45 100644
--- a/examples/solidHeatTransferScript/HeatConduction2DFin/HeatConduction2DFin.html
+++ b/examples/solidHeatTransferScript/HeatConduction2DFin/HeatConduction2DFin.html
@@ -21,29 +21,120 @@
     <script src="https://cdnjs.cloudflare.com/ajax/libs/plotly.js/2.27.0/plotly.min.js"></script>
 
     <!-- Link to the CSS files -->
-    <link href="https://feascript.com/FEAScript-website.css" rel="stylesheet" type="text/css" />
-    <link href="https://fonts.googleapis.com/css2?family=Roboto:wght@400;700&display=swap" rel="stylesheet" />
+    <link
+      href="https://feascript.com/FEAScript-website.css"
+      rel="stylesheet"
+      type="text/css"
+    />
+    <link
+      href="https://fonts.googleapis.com/css2?family=Roboto:wght@400;700&display=swap"
+      rel="stylesheet"
+    />
+
+    <!-- Custom layout styles -->
+    <style>
+      .container {
+        display: flex;
+        flex-wrap: wrap;
+      }
+
+      .column {
+        flex: 1 1 100%;
+        padding: 20px;
+        box-sizing: border-box;
+      }
+
+      @media (min-width: 768px) {
+        .column {
+          flex: 1 1 50%;
+        }
+      }
+      #solutionPlot {
+        width: 50%;
+        max-width: 100%;
+        min-height: 100px;
+        box-sizing: border-box;
+      }
+
+      #inputForm {
+        margin-top: 20px;
+      }
+      #resultOutput {
+        height: 70vh;
+        margin-top: 10px;
+        padding: 10px;
+        background: #f0f0f0;
+        overflow: auto;
+        border: 1px solid #ccc;
+        white-space: pre-wrap;
+      }
+    </style>
   </head>
 
   <body>
-    <h1>Heat Conduction in a Two-Dimensional Fin Example</h1>
-    <div id="solutionPlot"></div>
-
-    <p>
-      The mesh configuration and boundary conditions are defined directly within the JavaScript code in this
-      example. Please refresh the page to update the results. Detailed instructions for this example can be
-      found in the corresponding
-      <a href="https://feascript.com/tutorials/HeatConduction2DFin.html" target="_blank">FEAScript tutorial</a
-      >. If you need further assistance, you can visit the
-      <a href="https://feascript.com/" target="_blank">FEAScript website</a>.
-    </p>
-
-    <p>&#169; 2023-<span id="currentYear"></span> FEAScript</p>
-    <script>
-      document.getElementById("currentYear").innerHTML = new Date().getFullYear();
+    <div class="container">
+      <!-- Left Column: Simulation Content -->
+      <div class="column">
+        <h1>Heat Conduction in a Two-Dimensional Fin Example</h1>
+        <div id="solutionPlot"></div>
+
+        <p>
+          The mesh configuration and boundary conditions are defined directly
+          within the JavaScript code in this example. Please refresh the page to
+          update the results. Detailed instructions for this example can be
+          found in the corresponding
+          <a
+            href="https://feascript.com/tutorials/HeatConduction2DFin.html"
+            target="_blank"
+            >FEAScript tutorial</a
+          >. If you need further assistance, you can visit the
+          <a href="https://feascript.com/" target="_blank">FEAScript website</a
+          >.
+        </p>
+        <p>&#169; 2023-<span id="currentYear"></span> FEAScript</p>
+      </div>
+
+      <!-- Right Column: File Upload & Output -->
+      <div class="column">
+        <h2>Upload Input File</h2>
+        <form id="inputForm">
+          <label for="inputFile">Choose a file:</label>
+          <input
+            type="file"
+            id="inputFile"
+            name="inputFile"
+            accept=".msh"
+            required
+          />
+        </form>
+        <div id="resultOutput">No file uploaded.</div>
+      </div>
+    </div>
+
+    <script type="module">
+      import { importGmshQuadTri } from "../../../src/readers/gmshQuadReader.js";
+      document.getElementById("currentYear").innerHTML =
+        new Date().getFullYear();
+
+      document
+        .getElementById("inputFile")
+        .addEventListener("change", function (e) {
+          const file = e.target.files[0];
+          if (!file) return;
+          handleFileUpload(file);
+        });
+
+      async function handleFileUpload(file) {
+        const res = await importGmshQuadTri(file);
+        console.log(res);
+        document.getElementById("resultOutput").innerHTML = JSON.stringify(
+          res,
+          null,
+          2
+        );
+      }
     </script>
 
-    <!-- Import FEAScript library -->
     <script type="module">
       import {
         FEAScriptModel,
@@ -51,12 +142,14 @@ <h1>Heat Conduction in a Two-Dimensional Fin Example</h1>
         printVersion,
       } from "https://feascript.github.io/FEAScript-core/src/index.js";
 
-      window.addEventListener("DOMContentLoaded", () => {
+      import { FEAWorkerScript } from "../../../src/FEAWorkerScript.js";
+
+      window.addEventListener("DOMContentLoaded", async () => {
         // Print FEAScript version in the console
         printVersion();
 
         // Create a new FEAScript model
-        const model = new FEAScriptModel();
+        const model = new FEAWorkerScript();
 
         // Set solver configuration
         model.setSolverConfig("solidHeatTransferScript");
@@ -65,10 +158,10 @@ <h1>Heat Conduction in a Two-Dimensional Fin Example</h1>
         model.setMeshConfig({
           meshDimension: "2D",
           elementOrder: "quadratic",
-          numElementsX: 8,
-          numElementsY: 4,
-          maxX: 4,
-          maxY: 2,
+          numElementsX: 36,
+          numElementsY: 18,
+          maxX: 125,
+          maxY: 50,
         });
 
         // Define boundary conditions
@@ -81,14 +174,14 @@ <h1>Heat Conduction in a Two-Dimensional Fin Example</h1>
         model.setSolverMethod("lusolve");
 
         // Solve the problem and get the solution
-        const { solutionVector, nodesCoordinates } = model.solve();
+        const { solutionVector, nodesCoordinates } = await model.solve();
 
         // Plot the solution as a 2D contour plot
         plotSolution(
           solutionVector,
           nodesCoordinates,
           model.solverConfig,
-          model.meshConfig.meshDimension,
+          "2D",
           "contour",
           "solutionPlot"
         );
diff --git a/src/readers/gmshQuadReader.js b/src/readers/gmshQuadReader.js
index bca51b7..d2c7ab1 100644
--- a/src/readers/gmshQuadReader.js
+++ b/src/readers/gmshQuadReader.js
@@ -7,8 +7,9 @@ const importGmsh = async (file) => {
     gmshV: 0,
     ascii: false,
     fltBytes: 8,
-    totalNodesX:0,
-    totalNodesY:0
+    totalNodesX: 0,
+    totalNodesY: 0,
+    physicalPropMap: [],
   };
   let textre = await file.text();
   textre = textre
@@ -18,10 +19,17 @@ const importGmsh = async (file) => {
 
   let inNodesSections = false;
 
+  const totalElements = newFin.nodalNumbering.length;
+  const numMaterials = newFin.physicalPropMap.length;
+  const baseChunkSize = Math.floor(totalElements / numMaterials);
+  const remainder = totalElements % numMaterials;
+
   let inElementsSecitons = false;
 
   let inMeshSections = false;
 
+  let inPhysicalNames = false;
+
   let lineNumber = 0;
 
   for (let line of textre) {
@@ -34,6 +42,17 @@ const importGmsh = async (file) => {
       inMeshSections = false;
       continue;
     }
+    if (line == "$PhysicalNames") {
+      inPhysicalNames = true;
+      continue;
+    }
+    if (line == "$EndPhysicalNames") {
+      lineNumber = 0;
+      inPhysicalNames = false;
+
+      continue;
+    }
+
     if (line == "$Nodes") {
       inNodesSections = true;
       continue;
@@ -41,6 +60,7 @@ const importGmsh = async (file) => {
     if (line == "$EndNodes") {
       lineNumber = 0;
       inNodesSections = false;
+
       continue;
     }
     if (line == "$Elements") {
@@ -53,9 +73,10 @@ const importGmsh = async (file) => {
       inElementsSecitons = false;
       continue;
     }
+    let temp = line.split(" ");
+
     lineNumber = lineNumber + 1;
     if (inMeshSections) {
-      let temp = line.split(" ");
       let gmshVersion = parseFloat(temp[0]);
       let asciiOr = temp[1] === "0" ? true : false;
       let fltBytesOr = temp[2];
@@ -64,7 +85,23 @@ const importGmsh = async (file) => {
       newFin.ascii = asciiOr;
       newFin.fltBytes = fltBytesOr;
     }
-    let temp = line.split(" ");
+
+    if (inPhysicalNames) {
+      if (lineNumber > 1) {
+        let dimension = parseInt(temp[0]);
+        let tag = parseInt(temp[1]);
+        let name = temp[2].replace(/^"|"$/g, "");
+
+        newFin.physicalPropMap = [
+          ...newFin.physicalPropMap,
+          {
+            tag: tag,
+            dimension: dimension,
+            name: name,
+          },
+        ];
+      }
+    }
     if (inNodesSections) {
       if (temp.length === 3) {
         newFin.totalNodesX += 1;
@@ -90,64 +127,357 @@ const importGmsh = async (file) => {
     }
   }
 
-
   // const taichiRunner = new TaichiRunner();
   // await taichiRunner.init();
 
   // await taichiRunner.addToScope({ newFin, tempBoundary });
 
   // const useKernel = await taichiRunner.createKernel(() => {
-    const edgeCount = {};
-    for (let i = 0; i < newFin.nodalNumbering.length; i++) {
-      const element = newFin.nodalNumbering[i];
-      const edges = [
-        [element[0], element[1]],
-        [element[1], element[2]],
-        [element[2], element[3]],
-        [element[3], element[0]],
-      ];
-      edges.forEach((edge) => {
-        const key =
-          edge[0] < edge[1]
-            ? `${edge[0]}-${edge[1]}`
-            : `${edge[1]}-${edge[0]}`;
-        edgeCount[key] = (edgeCount[key] || 0) + 1;
-      });
-    }
-  
-    newFin.boundaryElements = []; 
-    for (let i = 0; i < newFin.nodalNumbering.length; i++) {
-      const element = newFin.nodalNumbering[i];
-      const edges = [
-        [element[0], element[1]],
-        [element[1], element[2]],
-        [element[2], element[3]],
-        [element[3], element[0]],
-      ];
-      for (let j = 0; j < edges.length; j++) {
-        const edge = edges[j];
-        const key =
-          edge[0] < edge[1]
-            ? `${edge[0]}-${edge[1]}`
-            : `${edge[1]}-${edge[0]}`;
-        if (edgeCount[key] === 1) {
-          if (!newFin.boundaryElements[i]) {
-            newFin.boundaryElements[i] = [];
-          }
-          newFin.boundaryElements[i].push([i, j]);
+  const edgeCount = {};
+  for (let i = 0; i < newFin.nodalNumbering.length; i++) {
+    const element = newFin.nodalNumbering[i];
+    const edges = [
+      [element[0], element[1]],
+      [element[1], element[2]],
+      [element[2], element[3]],
+      [element[3], element[0]],
+    ];
+    edges.forEach((edge) => {
+      const key =
+        edge[0] < edge[1] ? `${edge[0]}-${edge[1]}` : `${edge[1]}-${edge[0]}`;
+      edgeCount[key] = (edgeCount[key] || 0) + 1;
+    });
+  }
+
+  newFin.boundaryElements = [];
+  for (let i = 0; i < newFin.nodalNumbering.length; i++) {
+    const element = newFin.nodalNumbering[i];
+    const edges = [
+      [element[0], element[1]],
+      [element[1], element[2]],
+      [element[2], element[3]],
+      [element[3], element[0]],
+    ];
+    for (let j = 0; j < edges.length; j++) {
+      const edge = edges[j];
+      const key =
+        edge[0] < edge[1] ? `${edge[0]}-${edge[1]}` : `${edge[1]}-${edge[0]}`;
+      if (edgeCount[key] === 1) {
+        if (!newFin.boundaryElements[i]) {
+          newFin.boundaryElements[i] = [];
         }
+        newFin.boundaryElements[i].push([i, j]);
       }
     }
-  
-    newFin.boundaryElements = newFin.boundaryElements.filter((each) => each);
-  
+  }
+
+  newFin.boundaryElements = newFin.boundaryElements.filter((each) => each);
+
+  let tempElementArr = [];
+  let startIndex = 0;
+
+  for (let i = 0; i < numMaterials; i++) {
+    // Distribute one extra element for the first 'remainder' chunks
+    const currentChunkSize = baseChunkSize + (i < remainder ? 1 : 0);
+    tempElementArr.push(
+      newFin.nodalNumbering.slice(startIndex, startIndex + currentChunkSize)
+    );
+    startIndex += currentChunkSize;
+  }
+
+  newFin.nodalNumbering = tempElementArr;
+
   // });
 
   // if (useKernel) {
   //   await taichiRunner.runKernel(useKernel);
   // }
+  return newFin;
+};
+
+const importGmshQuadTri = async (file) => {
+  let newFin = {
+    nodesXCoordinates: [],
+    nodesYCoordinates: [],
+    nodalNumbering: {
+      quadElements: [],
+      triangleElements: [],
+    },
+    boundaryElements: [],
+    gmshV: 0,
+    ascii: false,
+    fltBytes: 8,
+    totalNodesX: 0,
+    totalNodesY: 0,
+    physicalPropMap: [],
+    elementTypes: {},
+  };
+
+  let textre = await file.text();
+  textre = textre
+    .split("\n")
+    .map((eachLine) => eachLine.trim())
+    .filter((eachLine) => eachLine != "" && eachLine != " ");
+
+  let inNodesSections = false;
+  let inElementsSecitons = false;
+  let inMeshSections = false;
+  let inPhysicalNames = false;
+  let inEntities = false;
+
+  let lineNumber = 0;
+  let entityBlocks = 0;
+  let totalNodes = 0;
+  let nodeBlocksProcessed = 0;
+  let currentNodeBlock = { dim: 0, tag: 0, parametric: 0, numNodes: 0 };
+  let nodeBlockDataStartLine = 0;
+  let collectingNodeTags = false;
+  let collectingNodeCoords = false;
+  let nodeTagsToProcess = [];
+
+  let entityBlocksElements = 0;
+  let totalElements = 0;
+  let elementBlocksProcessed = 0;
+  let currentElementBlock = { dim: 0, tag: 0, elementType: 0, numElements: 0 };
+  let elementBlockStartLine = 0;
+
+  for (let line of textre) {
+    if (line == "$MeshFormat") {
+      inMeshSections = true;
+      continue;
+    }
+    if (line == "$EndMeshFormat") {
+      lineNumber = 0;
+      inMeshSections = false;
+      continue;
+    }
+    if (line == "$PhysicalNames") {
+      inPhysicalNames = true;
+      continue;
+    }
+    if (line == "$EndPhysicalNames") {
+      lineNumber = 0;
+      inPhysicalNames = false;
+      continue;
+    }
+    if (line == "$Entities") {
+      inEntities = true;
+      continue;
+    }
+    if (line == "$EndEntities") {
+      lineNumber = 0;
+      inEntities = false;
+      continue;
+    }
+    if (line == "$Nodes") {
+      inNodesSections = true;
+      continue;
+    }
+    if (line == "$EndNodes") {
+      lineNumber = 0;
+      inNodesSections = false;
+      continue;
+    }
+    if (line == "$Elements") {
+      inElementsSecitons = true;
+      continue;
+    }
+    if (line == "$EndElements") {
+      lineNumber = 0;
+      inElementsSecitons = false;
+      continue;
+    }
+
+    let temp = line.split(" ").filter((item) => item !== "");
+
+    lineNumber = lineNumber + 1;
+
+    if (inMeshSections) {
+      if (lineNumber === 1) {
+        let gmshVersion = parseFloat(temp[0]);
+        let asciiOr = temp[1] === "0" ? true : false;
+        let fltBytesOr = temp[2];
+
+        newFin.gmshV = gmshVersion;
+        newFin.ascii = asciiOr;
+        newFin.fltBytes = fltBytesOr;
+      }
+    }
+
+    if (inPhysicalNames) {
+      if (lineNumber === 1) {
+        continue;
+      }
+
+      if (temp.length >= 3) {
+        let dimension = parseInt(temp[0]);
+        let tag = parseInt(temp[1]);
+
+        let name = temp.slice(2).join(" ");
+        name = name.replace(/^"|"$/g, "");
+
+        newFin.physicalPropMap.push({
+          tag: tag,
+          dimension: dimension,
+          name: name,
+        });
+      }
+    }
+
+    if (inNodesSections) {
+      if (lineNumber === 1) {
+        entityBlocks = parseInt(temp[0]);
+        totalNodes = parseInt(temp[1]);
+        newFin.nodesXCoordinates = new Array(totalNodes);
+        newFin.nodesYCoordinates = new Array(totalNodes);
+        continue;
+      }
+
+      if (nodeBlocksProcessed < entityBlocks) {
+        if (!collectingNodeTags && !collectingNodeCoords) {
+          currentNodeBlock.dim = parseInt(temp[0]);
+          currentNodeBlock.tag = parseInt(temp[1]);
+          currentNodeBlock.parametric = parseInt(temp[2]);
+          currentNodeBlock.numNodes = parseInt(temp[3]);
+
+          nodeTagsToProcess = [];
+          collectingNodeTags = true;
+          nodeBlockDataStartLine = lineNumber;
+          continue;
+        }
+
+        if (collectingNodeTags) {
+          for (let tag of temp) {
+            nodeTagsToProcess.push(parseInt(tag));
+          }
+
+          if (
+            lineNumber - nodeBlockDataStartLine >=
+            currentNodeBlock.numNodes
+          ) {
+            collectingNodeTags = false;
+            collectingNodeCoords = true;
+            continue;
+          }
+        }
+
+        if (collectingNodeCoords) {
+          const nodeIndex = nodeTagsToProcess.shift() - 1;
+
+          if (temp.length >= 2) {
+            newFin.nodesXCoordinates[nodeIndex] = parseFloat(temp[0]);
+            newFin.nodesYCoordinates[nodeIndex] = parseFloat(temp[1]);
+            newFin.totalNodesX++;
+            newFin.totalNodesY++;
+          }
+
+          if (nodeTagsToProcess.length === 0) {
+            collectingNodeCoords = false;
+            nodeBlocksProcessed++;
+          }
+        }
+      }
+    }
+
+    if (inElementsSecitons) {
+      if (lineNumber === 1) {
+        entityBlocksElements = parseInt(temp[0]);
+        totalElements = parseInt(temp[1]);
+        continue;
+      }
+
+      if (elementBlocksProcessed < entityBlocksElements) {
+        if (temp.length === 4) {
+          currentElementBlock.dim = parseInt(temp[0]);
+          currentElementBlock.tag = parseInt(temp[1]);
+          currentElementBlock.elementType = parseInt(temp[2]);
+          currentElementBlock.numElements = parseInt(temp[3]);
+
+          newFin.elementTypes[currentElementBlock.elementType] =
+            (newFin.elementTypes[currentElementBlock.elementType] || 0) +
+            currentElementBlock.numElements;
+
+          elementBlockStartLine = lineNumber;
+          continue;
+        }
+
+        if (
+          lineNumber > elementBlockStartLine &&
+          lineNumber <= elementBlockStartLine + currentElementBlock.numElements
+        ) {
+          const elementTag = parseInt(temp[0]);
+          const nodeIndices = temp.slice(1).map((idx) => parseInt(idx));
+
+          if (currentElementBlock.elementType === 3) {
+            newFin.nodalNumbering.quadElements.push(nodeIndices);
+          } else if (currentElementBlock.elementType === 2) {
+            newFin.nodalNumbering.triangleElements.push(nodeIndices);
+          }
+
+          if (
+            lineNumber ===
+            elementBlockStartLine + currentElementBlock.numElements
+          ) {
+            elementBlocksProcessed++;
+          }
+        }
+      }
+    }
+  }
+
+  processBoundaryElements(
+    newFin.nodalNumbering.triangleElements,
+    newFin.boundaryElements,
+    3
+  );
+  processBoundaryElements(
+    newFin.nodalNumbering.quadElements,
+    newFin.boundaryElements,
+    4
+  );
 
   return newFin;
 };
 
-export { importGmsh };
+function processBoundaryElements(elements, boundaryElements, numNodes) {
+  const edgeCount = {};
+
+  for (let i = 0; i < elements.length; i++) {
+    const element = elements[i];
+
+    const edges = [];
+    for (let j = 0; j < numNodes; j++) {
+      edges.push([element[j], element[(j + 1) % numNodes]]);
+    }
+
+    edges.forEach((edge) => {
+      const key =
+        edge[0] < edge[1] ? `${edge[0]}-${edge[1]}` : `${edge[1]}-${edge[0]}`;
+      edgeCount[key] = (edgeCount[key] || 0) + 1;
+    });
+  }
+  for (let i = 0; i < elements.length; i++) {
+    const element = elements[i];
+    const edges = [];
+
+    for (let j = 0; j < numNodes; j++) {
+      edges.push([element[j], element[(j + 1) % numNodes]]);
+    }
+
+    for (let j = 0; j < edges.length; j++) {
+      const edge = edges[j];
+      const key =
+        edge[0] < edge[1] ? `${edge[0]}-${edge[1]}` : `${edge[1]}-${edge[0]}`;
+
+      if (edgeCount[key] === 1) {
+        boundaryElements.push({
+          elementIndex: i,
+          localEdgeIndex: j,
+          elementType: numNodes === 3 ? "triangle" : "quad",
+        });
+      }
+    }
+  }
+}
+
+export { importGmsh, importGmshQuadTri };
diff --git a/src/visualization/plotSolutionScript.js b/src/visualization/plotSolutionScript.js
index dd873e5..67f1cf6 100644
--- a/src/visualization/plotSolutionScript.js
+++ b/src/visualization/plotSolutionScript.js
@@ -35,11 +35,20 @@ export function plotSolution(
     const numNodesY = new Set(nodesYCoordinates).size;
 
     // Reshape the nodesXCoordinates and nodesYCoordinates arrays to match the grid dimensions
-    let reshapedXCoordinates = math.reshape(Array.from(nodesXCoordinates), [numNodesX, numNodesY]);
-    let reshapedYCoordinates = math.reshape(Array.from(nodesYCoordinates), [numNodesX, numNodesY]);
+    let reshapedXCoordinates = math.reshape(Array.from(nodesXCoordinates), [
+      numNodesX,
+      numNodesY,
+    ]);
+    let reshapedYCoordinates = math.reshape(Array.from(nodesYCoordinates), [
+      numNodesX,
+      numNodesY,
+    ]);
 
     // Reshape the solution array to match the grid dimensions
-    let reshapedSolution = math.reshape(Array.from(solutionVector), [numNodesX, numNodesY]);
+    let reshapedSolution = math.reshape(Array.from(solutionVector), [
+      numNodesX,
+      numNodesY,
+    ]);
 
     // Transpose the reshapedSolution array to get column-wise data
     let transposedSolution = math.transpose(reshapedSolution);
@@ -104,7 +113,9 @@ export function plotSolution(
 
     // Set the layout for the contour plot
     let layout = {
-      title: `${plotType} plot${showMesh ? " with mesh" : ""} - ${solverConfig}`,
+      title: `${plotType} plot${
+        showMesh ? " with mesh" : ""
+      } - ${solverConfig}`,
       width: plotWidth,
       height: plotHeight,
       xaxis: { title: "x" },
@@ -116,6 +127,6 @@ export function plotSolution(
     if (showMesh) {
       plotData.push(meshData);
     }
-    Plotly.newPlot(plotDivId, plotData, layout);
+    Plotly.newPlot(plotDivId, plotData, layout, { responsive: true });
   }
 }
