feascript.esm.js

function e(e){let t=0;for(let n=0;n<e.length;n++)t+=e[n]*e[n];return t=Math.sqrt(t),t}let t="basic";function n(e){"basic"!==e&&"debug"!==e?(console.log("%c[WARN] Invalid log level: "+e+". Using basic instead.","color: #FFC107; font-weight: bold;"),t="basic"):(t=e,s(`Log level set to: ${e}`))}function o(e){"debug"===t&&console.log("%c[DEBUG] "+e,"color: #2196F3; font-weight: bold;")}function s(e){console.log("%c[INFO] "+e,"color: #4CAF50; font-weight: bold;")}function i(e){console.log("%c[ERROR] "+e,"color: #F44336; font-weight: bold;")}
/**
 * @license
 * Copyright 2019 Google LLC
 * SPDX-License-Identifier: Apache-2.0
 */
const r=Symbol("Comlink.proxy"),a=Symbol("Comlink.endpoint"),l=Symbol("Comlink.releaseProxy"),d=Symbol("Comlink.finalizer"),c=Symbol("Comlink.thrown"),u=e=>"object"==typeof e&&null!==e||"function"==typeof e,m=new Map([["proxy",{canHandle:e=>u(e)&&e[r],serialize(e){const{port1:t,port2:n}=new MessageChannel;return h(e,t),[n,[n]]},deserialize:e=>(e.start(),p(e))}],["throw",{canHandle:e=>u(e)&&c in e,serialize({value:e}){let t;return t=e instanceof Error?{isError:!0,value:{message:e.message,name:e.name,stack:e.stack}}:{isError:!1,value:e},[t,[]]},deserialize(e){if(e.isError)throw Object.assign(new Error(e.value.message),e.value);throw e.value}}]]);function h(e,t=globalThis,n=["*"]){t.addEventListener("message",(function o(s){if(!s||!s.data)return;if(!function(e,t){for(const n of e){if(t===n||"*"===n)return!0;if(n instanceof RegExp&&n.test(t))return!0}return!1}(n,s.origin))return void console.warn(`Invalid origin '${s.origin}' for comlink proxy`);const{id:i,type:a,path:l}=Object.assign({path:[]},s.data),u=(s.data.argumentList||[]).map(F);let m;try{const t=l.slice(0,-1).reduce(((e,t)=>e[t]),e),n=l.reduce(((e,t)=>e[t]),e);switch(a){case"GET":m=n;break;case"SET":t[l.slice(-1)[0]]=F(s.data.value),m=!0;break;case"APPLY":m=n.apply(t,u);break;case"CONSTRUCT":m=function(e){return Object.assign(e,{[r]:!0})}(new n(...u));break;case"ENDPOINT":{const{port1:t,port2:n}=new MessageChannel;h(e,n),m=function(e,t){return D.set(e,t),e}(t,[t])}break;case"RELEASE":m=void 0;break;default:return}}catch(e){m={value:e,[c]:0}}Promise.resolve(m).catch((e=>({value:e,[c]:0}))).then((n=>{const[s,r]=M(n);t.postMessage(Object.assign(Object.assign({},s),{id:i}),r),"RELEASE"===a&&(t.removeEventListener("message",o),f(t),d in e&&"function"==typeof e[d]&&e[d]())})).catch((e=>{const[n,o]=M({value:new TypeError("Unserializable return value"),[c]:0});t.postMessage(Object.assign(Object.assign({},n),{id:i}),o)}))})),t.start&&t.start()}function f(e){(function(e){return"MessagePort"===e.constructor.name})(e)&&e.close()}function p(e,t){const n=new Map;return e.addEventListener("message",(function(e){const{data:t}=e;if(!t||!t.id)return;const o=n.get(t.id);if(o)try{o(t)}finally{n.delete(t.id)}})),E(e,n,[],t)}function g(e){if(e)throw new Error("Proxy has been released and is not useable")}function b(e){return $(e,new Map,{type:"RELEASE"}).then((()=>{f(e)}))}const y=new WeakMap,v="FinalizationRegistry"in globalThis&&new FinalizationRegistry((e=>{const t=(y.get(e)||0)-1;y.set(e,t),0===t&&b(e)}));function E(e,t,n=[],o=function(){}){let s=!1;const i=new Proxy(o,{get(o,r){if(g(s),r===l)return()=>{!function(e){v&&v.unregister(e)}(i),b(e),t.clear(),s=!0};if("then"===r){if(0===n.length)return{then:()=>i};const o=$(e,t,{type:"GET",path:n.map((e=>e.toString()))}).then(F);return o.then.bind(o)}return E(e,t,[...n,r])},set(o,i,r){g(s);const[a,l]=M(r);return $(e,t,{type:"SET",path:[...n,i].map((e=>e.toString())),value:a},l).then(F)},apply(o,i,r){g(s);const l=n[n.length-1];if(l===a)return $(e,t,{type:"ENDPOINT"}).then(F);if("bind"===l)return E(e,t,n.slice(0,-1));const[d,c]=C(r);return $(e,t,{type:"APPLY",path:n.map((e=>e.toString())),argumentList:d},c).then(F)},construct(o,i){g(s);const[r,a]=C(i);return $(e,t,{type:"CONSTRUCT",path:n.map((e=>e.toString())),argumentList:r},a).then(F)}});return function(e,t){const n=(y.get(t)||0)+1;y.set(t,n),v&&v.register(e,t,e)}(i,e),i}function C(e){const t=e.map(M);return[t.map((e=>e[0])),(n=t.map((e=>e[1])),Array.prototype.concat.apply([],n))];var n}const D=new WeakMap;function M(e){for(const[t,n]of m)if(n.canHandle(e)){const[o,s]=n.serialize(e);return[{type:"HANDLER",name:t,value:o},s]}return[{type:"RAW",value:e},D.get(e)||[]]}function F(e){switch(e.type){case"HANDLER":return m.get(e.name).deserialize(e.value);case"RAW":return e.value}}function $(e,t,n,o){return new Promise((s=>{const i=new Array(4).fill(0).map((()=>Math.floor(Math.random()*Number.MAX_SAFE_INTEGER).toString(16))).join("-");t.set(i,s),e.start&&e.start(),e.postMessage(Object.assign({id:i},n),o)}))}function w(e,t,n,r={}){const{maxIterations:a=1e4,tolerance:l=1e-4}=r;let d=[],c=!0,u=0;if(s(`Solving system using ${e}...`),console.time("systemSolving"),"lusolve"===e){const e=math.sparse(t),o=math.slu(e,1,1);let s=math.lusolve(o,n);d=math.squeeze(s).valueOf()}else if("jacobi"===e){const e=function(e,t,n,o={}){const{maxIterations:s,tolerance:i}=o,r=e.length;let a=[...n],l=new Array(r);for(let n=0;n<s;n++){for(let n=0;n<r;n++){let o=0;for(let t=0;t<r;t++)n!==t&&(o+=e[n][t]*a[t]);l[n]=(t[n]-o)/e[n][n]}let o=0;for(let e=0;e<r;e++)o=Math.max(o,Math.abs(l[e]-a[e]));if(a=[...l],o<i)return{solutionVector:a,iterations:n+1,converged:!0}}return{solutionVector:a,iterations:s,converged:!1}}(t,n,new Array(n.length).fill(0),{maxIterations:a,tolerance:l});e.converged?o(`Jacobi method converged in ${e.iterations} iterations`):i(`Jacobi method did not converge after ${e.iterations} iterations`),d=e.solutionVector,c=e.converged,u=e.iterations}else i(`Unknown solver method: ${e}`);return console.timeEnd("systemSolving"),s("System solved successfully"),{solutionVector:d,converged:c,iterations:u}}async function A(e,t,n,r={}){const{maxIterations:a=1e4,tolerance:l=1e-4}=r;s(`Solving system using ${e}...`),console.time("systemSolving");const d=Array.isArray(t)?t:t?.toArray?.()??t,c=Array.isArray(n)?n:n?.toArray?.()??n;let u,m=null,h=null,f=[],g=!0;if("jacobi-gpu"===e){m=await async function(){const e=new Worker(new URL("../workers/webgpuWorkerScript.js",import.meta.url),{type:"module"}),t=p(e);return await t.initialize(),{computeEngine:t,worker:e}}(),h=m.computeEngine;const e=new Array(c.length).fill(0);let t;t=await h.webgpuJacobiSolver(d,c,e,{maxIterations:a,tolerance:l}),f=t.solutionVector,g=t.converged,u=t.iterations,g?o(`Jacobi method converged in ${u} iterations`):i(`Jacobi method did not converge after ${u} iterations`)}else i(`Unknown solver method: ${e}`);return console.timeEnd("systemSolving"),s(`System solved successfully (${e})`),m&&(await(h?.destroy?.().catch((()=>{}))),m.worker.terminate()),{solutionVector:f,converged:g,iterations:u}}class x{constructor({meshDimension:e,elementOrder:t}){this.meshDimension=e,this.elementOrder=t}getBasisFunctions(e,t=null){let n=[],o=[],s=[];if("1D"===this.meshDimension)"linear"===this.elementOrder?(n[0]=1-e,n[1]=e,o[0]=-1,o[1]=1):"quadratic"===this.elementOrder&&(n[0]=1-3*e+2*e**2,n[1]=4*e-4*e**2,n[2]=2*e**2-e,o[0]=4*e-3,o[1]=4-8*e,o[2]=4*e-1);else if("2D"===this.meshDimension){if(null===t)return void i("Eta coordinate is required for 2D elements");if("linear"===this.elementOrder){function r(e){return 1-e}n[0]=r(e)*r(t),n[1]=r(e)*t,n[2]=e*r(t),n[3]=e*t,o[0]=-1*r(t),o[1]=-1*t,o[2]=1*r(t),o[3]=1*t,s[0]=-1*r(e),s[1]=1*r(e),s[2]=-1*e,s[3]=1*e}else if("quadratic"===this.elementOrder){function a(e){return 2*e**2-3*e+1}function l(e){return-4*e**2+4*e}function d(e){return 2*e**2-e}function c(e){return 4*e-3}function u(e){return-8*e+4}function m(e){return 4*e-1}n[0]=a(e)*a(t),n[1]=a(e)*l(t),n[2]=a(e)*d(t),n[3]=l(e)*a(t),n[4]=l(e)*l(t),n[5]=l(e)*d(t),n[6]=d(e)*a(t),n[7]=d(e)*l(t),n[8]=d(e)*d(t),o[0]=c(e)*a(t),o[1]=c(e)*l(t),o[2]=c(e)*d(t),o[3]=u(e)*a(t),o[4]=u(e)*l(t),o[5]=u(e)*d(t),o[6]=m(e)*a(t),o[7]=m(e)*l(t),o[8]=m(e)*d(t),s[0]=a(e)*c(t),s[1]=a(e)*u(t),s[2]=a(e)*m(t),s[3]=l(e)*c(t),s[4]=l(e)*u(t),s[5]=l(e)*m(t),s[6]=d(e)*c(t),s[7]=d(e)*u(t),s[8]=d(e)*m(t)}}return{basisFunction:n,basisFunctionDerivKsi:o,basisFunctionDerivEta:s}}}class N{constructor({numElementsX:e=null,maxX:t=null,numElementsY:n=null,maxY:o=null,meshDimension:i=null,elementOrder:r="linear",parsedMesh:a=null}){this.numElementsX=e,this.numElementsY=n,this.maxX=t,this.maxY=o,this.meshDimension=i,this.elementOrder=r,this.parsedMesh=a,this.boundaryElementsProcessed=!1,this.parsedMesh&&(s("Using pre-parsed mesh from gmshReader data for mesh generation."),this.parseMeshFromGmsh())}parseMeshFromGmsh(){if(this.parsedMesh.nodalNumbering||i("No valid nodal numbering found in the parsed mesh."),Array.isArray(this.parsedMesh.nodalNumbering))return this.boundaryElementsProcessed=!0,this.parsedMesh.boundaryElementsProcessed=!0,this.parsedMesh;if("object"==typeof this.parsedMesh.nodalNumbering&&!Array.isArray(this.parsedMesh.nodalNumbering)){const e=this.parsedMesh.nodalNumbering.quadElements||[];if(this.parsedMesh.nodalNumbering.triangleElements,o("Initial parsed mesh nodal numbering from Gmsh format: "+JSON.stringify(this.parsedMesh.nodalNumbering)),this.parsedMesh.elementTypes[3]||this.parsedMesh.elementTypes[10]){const t=[];for(let n=0;n<e.length;n++){const o=e[n],s=new Array(o.length);4===o.length?(s[0]=o[0],s[1]=o[3],s[2]=o[1],s[3]=o[2]):9===o.length&&(s[0]=o[0],s[1]=o[7],s[2]=o[3],s[3]=o[4],s[4]=o[8],s[5]=o[6],s[6]=o[1],s[7]=o[5],s[8]=o[2]),t.push(s)}this.parsedMesh.nodalNumbering=t}else this.parsedMesh.elementTypes[2]&&i("Element type is neither triangle nor quad; mapping for this type is not implemented yet.");if(o("Nodal numbering after mapping from Gmsh to FEAScript format: "+JSON.stringify(this.parsedMesh.nodalNumbering)),this.parsedMesh.physicalPropMap&&this.parsedMesh.boundaryElements){if(Array.isArray(this.parsedMesh.boundaryElements)&&this.parsedMesh.boundaryElements.length>0&&void 0===this.parsedMesh.boundaryElements[0]){const e=[];for(let t=1;t<this.parsedMesh.boundaryElements.length;t++)this.parsedMesh.boundaryElements[t]&&e.push(this.parsedMesh.boundaryElements[t]);this.parsedMesh.boundaryElements=e}if(this.parsedMesh.boundaryNodePairs&&!this.parsedMesh.boundaryElementsProcessed&&(this.parsedMesh.boundaryElements=[],this.parsedMesh.physicalPropMap.forEach((e=>{if(1===e.dimension){const t=this.parsedMesh.boundaryNodePairs[e.tag]||[];t.length>0&&(this.parsedMesh.boundaryElements[e.tag]||(this.parsedMesh.boundaryElements[e.tag]=[]),t.forEach((t=>{const n=t[0],s=t[1];o(`Processing boundary node pair: [${n}, ${s}] for boundary ${e.tag} (${e.name||"unnamed"})`);let r=!1;for(let t=0;t<this.parsedMesh.nodalNumbering.length;t++){const i=this.parsedMesh.nodalNumbering[t];if(4===i.length){if(i.includes(n)&&i.includes(s)){let a;const l=i.indexOf(n),d=i.indexOf(s);o(`  Found element ${t} containing boundary nodes. Element nodes: [${i.join(", ")}]`),o(`  Node ${n} is at index ${l}, Node ${s} is at index ${d} in the element`),0===l&&2===d||2===l&&0===d?(a=0,o(`  These nodes form the BOTTOM side (${a}) of element ${t}`)):0===l&&1===d||1===l&&0===d?(a=1,o(`  These nodes form the LEFT side (${a}) of element ${t}`)):1===l&&3===d||3===l&&1===d?(a=2,o(`  These nodes form the TOP side (${a}) of element ${t}`)):(2===l&&3===d||3===l&&2===d)&&(a=3,o(`  These nodes form the RIGHT side (${a}) of element ${t}`)),this.parsedMesh.boundaryElements[e.tag].push([t,a]),o(`  Added element-side pair [${t}, ${a}] to boundary tag ${e.tag}`),r=!0;break}}else if(9===i.length&&i.includes(n)&&i.includes(s)){let a;const l=i.indexOf(n),d=i.indexOf(s);o(`  Found element ${t} containing boundary nodes. Element nodes: [${i.join(", ")}]`),o(`  Node ${n} is at index ${l}, Node ${s} is at index ${d} in the element`),0===l&&6===d||6===l&&0===d||0===l&&3===d||3===l&&0===d||3===l&&6===d||6===l&&3===d?(a=0,o(`  These nodes form the BOTTOM side (${a}) of element ${t}`)):0===l&&2===d||2===l&&0===d||0===l&&1===d||1===l&&0===d||1===l&&2===d||2===l&&1===d?(a=1,o(`  These nodes form the LEFT side (${a}) of element ${t}`)):2===l&&8===d||8===l&&2===d||2===l&&5===d||5===l&&2===d||5===l&&8===d||8===l&&5===d?(a=2,o(`  These nodes form the TOP side (${a}) of element ${t}`)):(6===l&&8===d||8===l&&6===d||6===l&&7===d||7===l&&6===d||7===l&&8===d||8===l&&7===d)&&(a=3,o(`  These nodes form the RIGHT side (${a}) of element ${t}`)),this.parsedMesh.boundaryElements[e.tag].push([t,a]),o(`  Added element-side pair [${t}, ${a}] to boundary tag ${e.tag}`),r=!0;break}}r||i(`Could not find element containing boundary nodes ${n} and ${s}. Boundary may be incomplete.`)})))}})),this.boundaryElementsProcessed=!0,this.parsedMesh.boundaryElements.length>0&&void 0===this.parsedMesh.boundaryElements[0])){const e=[];for(let t=1;t<this.parsedMesh.boundaryElements.length;t++)this.parsedMesh.boundaryElements[t]&&e.push(this.parsedMesh.boundaryElements[t]);this.parsedMesh.boundaryElements=e}}}return this.parsedMesh}}class k extends N{constructor({numElementsX:e=null,maxX:t=null,elementOrder:n="linear",parsedMesh:o=null}){super({numElementsX:e,maxX:t,numElementsY:1,maxY:0,meshDimension:"1D",elementOrder:n,parsedMesh:o}),null!==this.numElementsX&&null!==this.maxX||i("numElementsX and maxX are required parameters when generating a 1D mesh from geometry")}generateMesh(){let e=[];let t,n;if("linear"===this.elementOrder){t=this.numElementsX+1,n=(this.maxX-0)/this.numElementsX,e[0]=0;for(let o=1;o<t;o++)e[o]=e[o-1]+n}else if("quadratic"===this.elementOrder){t=2*this.numElementsX+1,n=(this.maxX-0)/this.numElementsX,e[0]=0;for(let o=1;o<t;o++)e[o]=e[o-1]+n/2}const s=this.generateNodalNumbering1D(this.numElementsX,t,this.elementOrder),i=this.findBoundaryElements();return o("Generated node X coordinates: "+JSON.stringify(e)),{nodesXCoordinates:e,totalNodesX:t,nodalNumbering:s,boundaryElements:i}}generateNodalNumbering1D(e,t,n){let o=[];if("linear"===n)for(let t=0;t<e;t++){o[t]=[];for(let e=1;e<=2;e++)o[t][e-1]=t+e}else if("quadratic"===n){let t=0;for(let n=0;n<e;n++){o[n]=[];for(let e=1;e<=3;e++)o[n][e-1]=n+e+t;t+=1}}return o}findBoundaryElements(){const e=[];for(let t=0;t<2;t++)e.push([]);return e[0].push([0,0]),e[1].push([this.numElementsX-1,1]),o("Identified boundary elements by side: "+JSON.stringify(e)),this.boundaryElementsProcessed=!0,e}}class P extends N{constructor({numElementsX:e=null,maxX:t=null,numElementsY:n=null,maxY:o=null,elementOrder:s="linear",parsedMesh:r=null,angleLeft:a=90,angleRight:l=90}){super({numElementsX:e,maxX:t,numElementsY:n,maxY:o,meshDimension:"2D",elementOrder:s,parsedMesh:r}),this.angleLeft=a,this.angleRight=l,r||null!==this.numElementsX&&null!==this.maxX&&null!==this.numElementsY&&null!==this.maxY||i("numElementsX, maxX, numElementsY, and maxY are required parameters when generating a 2D mesh from geometry")}generateMesh(){let e=[],t=[];let n,s,i,r;if("linear"===this.elementOrder){n=this.numElementsX+1,s=this.numElementsY+1,i=(this.maxX-0)/this.numElementsX,r=(this.maxY-0)/this.numElementsY,e[0]=0,t[0]=0;for(let n=1;n<s;n++)e[n]=e[0],t[n]=t[0]+n*r;for(let o=1;o<n;o++){const n=o*s;e[n]=e[0]+o*i,t[n]=t[0];for(let o=1;o<s;o++)e[n+o]=e[n],t[n+o]=t[n]+o*r}}else if("quadratic"===this.elementOrder){n=2*this.numElementsX+1,s=2*this.numElementsY+1,i=(this.maxX-0)/this.numElementsX,r=(this.maxY-0)/this.numElementsY,e[0]=0,t[0]=0;for(let n=1;n<s;n++)e[n]=e[0],t[n]=t[0]+n*r/2;for(let o=1;o<n;o++){const n=o*s;e[n]=e[0]+o*i/2,t[n]=t[0];for(let o=1;o<s;o++)e[n+o]=e[n],t[n+o]=t[n]+o*r/2}}if(90!==this.angleLeft||90!==this.angleRight){const o=Math.PI/180,i=Math.tan(this.angleLeft*o),r=Math.tan(this.angleRight*o),a=1e-12;for(let o=0;o<s;o++){const l=t[o],d=Math.abs(i)<a?0:l/i,c=Math.abs(r)<a?this.maxX:this.maxX-l/r;for(let t=0;t<n;t++){const i=1===n?0:t/(n-1);e[t*s+o]=d+(c-d)*i}}}const a=this.generateNodalNumbering2D(this.numElementsX,this.numElementsY,s,this.elementOrder),l=this.findBoundaryElements();return o("Generated node X coordinates: "+JSON.stringify(e)),o("Generated node Y coordinates: "+JSON.stringify(t)),{nodesXCoordinates:e,nodesYCoordinates:t,totalNodesX:n,totalNodesY:s,nodalNumbering:a,boundaryElements:l}}generateNodalNumbering2D(e,t,n,o){let s=0,i=[];if("linear"===o){let n=0,o=2;for(let s=0;s<e*t;s++)n+=1,i[s]=[],i[s][0]=s+o-1,i[s][1]=s+o,i[s][2]=s+o+t,i[s][3]=s+o+t+1,n===t&&(o+=1,n=0)}else if("quadratic"===o)for(let o=1;o<=e;o++)for(let e=1;e<=t;e++){i[s]=[];for(let t=1;t<=3;t++){let r=3*t-2;i[s][r-1]=n*(2*o+t-3)+2*e-1,i[s][r]=i[s][r-1]+1,i[s][r+1]=i[s][r-1]+2}s+=1}return i}findBoundaryElements(){const e=[];for(let t=0;t<4;t++)e.push([]);for(let t=0;t<this.numElementsX;t++)for(let n=0;n<this.numElementsY;n++){const o=t*this.numElementsY+n;0===n&&e[0].push([o,0]),0===t&&e[1].push([o,1]),n===this.numElementsY-1&&e[2].push([o,2]),t===this.numElementsX-1&&e[3].push([o,3])}return o("Identified boundary elements by side: "+JSON.stringify(e)),this.boundaryElementsProcessed=!0,e}}class O{constructor({meshDimension:e,elementOrder:t}){this.meshDimension=e,this.elementOrder=t}getGaussPointsAndWeights(){let e=[],t=[];return"linear"===this.elementOrder?(e[0]=.5,t[0]=1):"quadratic"===this.elementOrder&&(e[0]=(1-Math.sqrt(.6))/2,e[1]=.5,e[2]=(1+Math.sqrt(.6))/2,t[0]=5/18,t[1]=8/18,t[2]=5/18),{gaussPoints:e,gaussWeights:t}}}function V(e){const{meshDimension:t,numElementsX:n,numElementsY:s,maxX:r,maxY:a,elementOrder:l,parsedMesh:d,angleLeft:c,angleRight:u}=e;let m;"1D"===t?m=new k({numElementsX:n,maxX:r,elementOrder:l,parsedMesh:d}):"2D"===t?m=new P({numElementsX:n,maxX:r,numElementsY:s,maxY:a,elementOrder:l,parsedMesh:d,angleLeft:c,angleRight:u}):i("Mesh dimension must be either '1D' or '2D'");const h=m.boundaryElementsProcessed?m.parsedMesh:m.generateMesh();let f=h.nodesXCoordinates,p=h.nodesYCoordinates,g=h.totalNodesX,b=h.totalNodesY,y=h.nodalNumbering,v=h.boundaryElements;let E,C;return null!=d?(E=y.length,C=f.length,o(`Using parsed mesh with ${E} elements and ${C} nodes`)):(E=n*("2D"===t?s:1),C=g*("2D"===t?b:1),o(`Using mesh generated from geometry with ${E} elements and ${C} nodes`)),{nodesXCoordinates:f,nodesYCoordinates:p,totalNodesX:g,totalNodesY:b,nop:y,boundaryElements:v,totalElements:E,totalNodes:C,meshDimension:t,elementOrder:l}}function X(e){const{totalNodes:t,nop:n,meshDimension:o,elementOrder:s}=e;let i=[],r=[];for(let e=0;e<t;e++){i[e]=0,r.push([]);for(let n=0;n<t;n++)r[e][n]=0}const a=new x({meshDimension:o,elementOrder:s});let l=new O({meshDimension:o,elementOrder:s}).getGaussPointsAndWeights();return{residualVector:i,jacobianMatrix:r,localToGlobalMap:[],basisFunctions:a,gaussPoints:l.gaussPoints,gaussWeights:l.gaussWeights,nodesPerElement:n[0].length}}function I(e){const{basisFunction:t,basisFunctionDerivKsi:n,nodesXCoordinates:o,localToGlobalMap:s,nodesPerElement:i}=e;let r=0,a=0;for(let e=0;e<i;e++)r+=o[s[e]]*t[e],a+=o[s[e]]*n[e];let l=a,d=[];for(let e=0;e<i;e++)d[e]=n[e]/l;return{xCoordinates:r,detJacobian:l,basisFunctionDerivX:d}}function S(e){const{basisFunction:t,basisFunctionDerivKsi:n,basisFunctionDerivEta:o,nodesXCoordinates:s,nodesYCoordinates:i,localToGlobalMap:r,nodesPerElement:a}=e;let l=0,d=0,c=0,u=0,m=0,h=0;for(let e=0;e<a;e++)l+=s[r[e]]*t[e],d+=i[r[e]]*t[e],c+=s[r[e]]*n[e],u+=s[r[e]]*o[e],m+=i[r[e]]*n[e],h+=i[r[e]]*o[e];let f=c*h-u*m,p=[],g=[];for(let e=0;e<a;e++)p[e]=(h*n[e]-m*o[e])/f,g[e]=(c*o[e]-u*n[e])/f;return{xCoordinates:l,yCoordinates:d,detJacobian:f,basisFunctionDerivX:p,basisFunctionDerivY:g}}function T(e,t,n){const[o,s,i]=n,r=(s[1]-i[1])*(o[0]-i[0])+(i[0]-s[0])*(o[1]-i[1]),a=((s[1]-i[1])*(e-i[0])+(i[0]-s[0])*(t-i[1]))/r,l=((i[1]-o[1])*(e-i[0])+(o[0]-i[0])*(t-i[1]))/r;return{inside:a>=-1e-12&&l>=-1e-12&&1-a-l>=-1e-12,ksi:a,eta:l}}function R(e,t,n){const[o,s]=function(e){const[t,n,o,s]=e;return[[t,n,s],[t,o,s]]}(n),i=T(e,t,o),r=T(e,t,s),a=i.inside||r.inside;let l=0,d=0;if(a){const[o,s,i,r]=n,a=(n,o)=>Math.abs((o[0]-n[0])*(n[1]-t)-(n[0]-e)*(o[1]-n[1]))/Math.sqrt((o[0]-n[0])**2+(o[1]-n[1])**2),c=a(o,s),u=a(i,r),m=a(o,i);l=c/(c+u),d=m/(m+a(s,r))}return{inside:a,ksi:l,eta:d}}function Y(e){const{nop:t,nodesXCoordinates:n}=e,o=n.length,s=t[0].length,i=Array.from({length:o},(()=>[])),r=Array(o).fill(0);for(let e=0;e<t.length;e++)for(let n=0;n<s;n++){const o=t[e][n]-1;r[o]=r[o]+1,i[o].push(e)}return{nodeNeighbors:i,neighborCount:r}}function B(e){let t,n=[],o=[],s=0;const{nodesXCoordinates:i,nodesYCoordinates:r,nop:a,boundaryElements:l,meshDimension:d,elementOrder:c}=e;"1D"===d?("linear"===c||"quadratic"===c)&&(t={0:[0],1:[1]}):"2D"===d&&("linear"===c?t={0:[0,2],1:[0,1],2:[1,3],3:[2,3]}:"quadratic"===c&&(t={0:[0,3,6],1:[0,1,2],2:[2,5,8],3:[6,7,8]}));for(let e=0;e<l.length;e++)for(let d=0;d<l[e].length;d++){n[s]=l[e][d],s++;const[c,u]=l[e][d];let m=t[u],h=[],f=[];for(let e=0;e<m.length;e++){const t=a[c][m[e]]-1;h.push(i[t]),f.push(r[t])}for(let e=0;e<h.length-1;e++)o.push([[h[e],f[e]],[h[e+1],f[e+1]]])}return o}class j{constructor(e,t,n,o,s){this.boundaryConditions=e,this.boundaryElements=t,this.nop=n,this.meshDimension=o,this.elementOrder=s}imposeConstantTempBoundaryConditions(e,t){"1D"===this.meshDimension?Object.keys(this.boundaryConditions).forEach((n=>{if("constantTemp"===this.boundaryConditions[n][0]){const s=this.boundaryConditions[n][1];o(`Boundary ${n}: Applying constant temperature of ${s} K (Dirichlet condition)`),this.boundaryElements[n].forEach((([n,i])=>{if("linear"===this.elementOrder){({0:[0],1:[1]})[i].forEach((i=>{const r=this.nop[n][i]-1;o(`  - Applied constant temperature to node ${r+1} (element ${n+1}, local node ${i+1})`),e[r]=s;for(let n=0;n<e.length;n++)t[r][n]=0;t[r][r]=1}))}else if("quadratic"===this.elementOrder){({0:[0],1:[2]})[i].forEach((i=>{const r=this.nop[n][i]-1;o(`  - Applied constant temperature to node ${r+1} (element ${n+1}, local node ${i+1})`),e[r]=s;for(let n=0;n<e.length;n++)t[r][n]=0;t[r][r]=1}))}}))}})):"2D"===this.meshDimension&&Object.keys(this.boundaryConditions).forEach((n=>{if("constantTemp"===this.boundaryConditions[n][0]){const s=this.boundaryConditions[n][1];o(`Boundary ${n}: Applying constant temperature of ${s} K (Dirichlet condition)`),this.boundaryElements[n].forEach((([n,i])=>{if("linear"===this.elementOrder){({0:[0,2],1:[0,1],2:[1,3],3:[2,3]})[i].forEach((i=>{const r=this.nop[n][i]-1;o(`  - Applied constant temperature to node ${r+1} (element ${n+1}, local node ${i+1})`),e[r]=s;for(let n=0;n<e.length;n++)t[r][n]=0;t[r][r]=1}))}else if("quadratic"===this.elementOrder){({0:[0,3,6],1:[0,1,2],2:[2,5,8],3:[6,7,8]})[i].forEach((i=>{const r=this.nop[n][i]-1;o(`  - Applied constant temperature to node ${r+1} (element ${n+1}, local node ${i+1})`),e[r]=s;for(let n=0;n<e.length;n++)t[r][n]=0;t[r][r]=1}))}}))}}))}imposeConstantTempBoundaryConditionsFront(e,t){"1D"===this.meshDimension?Object.keys(this.boundaryConditions).forEach((n=>{if("constantTemp"===this.boundaryConditions[n][0]){const s=this.boundaryConditions[n][1];o(`Boundary ${n}: Applying constant temperature of ${s} K (Dirichlet condition)`),this.boundaryElements[n].forEach((([n,i])=>{if("linear"===this.elementOrder){({0:[0],1:[1]})[i].forEach((i=>{const r=this.nop[n][i]-1;o(`  - Applied constant temperature to node ${r+1} (element ${n+1}, local node ${i+1})`),e[r]=1,t[r]=s}))}else if("quadratic"===this.elementOrder){({0:[0],2:[2]})[i].forEach((i=>{const r=this.nop[n][i]-1;o(`  - Applied constant temperature to node ${r+1} (element ${n+1}, local node ${i+1})`),e[r]=1,t[r]=s}))}}))}})):"2D"===this.meshDimension&&Object.keys(this.boundaryConditions).forEach((n=>{if("constantTemp"===this.boundaryConditions[n][0]){const s=this.boundaryConditions[n][1];o(`Boundary ${n}: Applying constant temperature of ${s} K (Dirichlet condition)`),this.boundaryElements[n].forEach((([n,i])=>{if("linear"===this.elementOrder){({0:[0,2],1:[0,1],2:[1,3],3:[2,3]})[i].forEach((i=>{const r=this.nop[n][i]-1;o(`  - Applied constant temperature to node ${r+1} (element ${n+1}, local node ${i+1})`),e[r]=1,t[r]=s}))}else if("quadratic"===this.elementOrder){({0:[0,3,6],1:[0,1,2],2:[2,5,8],3:[6,7,8]})[i].forEach((i=>{const r=this.nop[n][i]-1;o(`  - Applied constant temperature to node ${r+1} (element ${n+1}, local node ${i+1})`),e[r]=1,t[r]=s}))}}))}}))}imposeConvectionBoundaryConditions(e,t,n,s,i,r,a){let l=[],d=[];Object.keys(this.boundaryConditions).forEach((e=>{const t=this.boundaryConditions[e];"convection"===t[0]&&(l[e]=t[1],d[e]=t[2])})),"1D"===this.meshDimension?Object.keys(this.boundaryConditions).forEach((n=>{if("convection"===this.boundaryConditions[n][0]){const s=l[n],i=d[n];o(`Boundary ${n}: Applying convection with heat transfer coefficient h=${s} W/(m²·K) and external temperature T∞=${i} K`),this.boundaryElements[n].forEach((([n,r])=>{let a;"linear"===this.elementOrder?a=0===r?0:1:"quadratic"===this.elementOrder&&(a=0===r?0:2);const l=this.nop[n][a]-1;o(`  - Applied convection boundary condition to node ${l+1} (element ${n+1}, local node ${a+1})`),e[l]+=-s*i,t[l][l]+=s}))}})):"2D"===this.meshDimension&&Object.keys(this.boundaryConditions).forEach((c=>{if("convection"===this.boundaryConditions[c][0]){const u=l[c],m=d[c];o(`Boundary ${c}: Applying convection with heat transfer coefficient h=${u} W/(m²·K) and external temperature T∞=${m} K`),this.boundaryElements[c].forEach((([l,d])=>{if("linear"===this.elementOrder){let c,h,f,p,g;0===d?(c=n[0],h=0,f=0,p=3,g=2):1===d?(c=0,h=n[0],f=0,p=2,g=1):2===d?(c=n[0],h=1,f=1,p=4,g=2):3===d&&(c=1,h=n[0],f=2,p=4,g=1);let b=a.getBasisFunctions(c,h),y=b.basisFunction,v=b.basisFunctionDerivKsi,E=b.basisFunctionDerivEta,C=0,D=0,M=0,F=0;const $=this.nop[l].length;for(let e=0;e<$;e++){const t=this.nop[l][e]-1;0===d||2===d?(C+=i[t]*v[e],D+=r[t]*v[e]):1!==d&&3!==d||(M+=i[t]*E[e],F+=r[t]*E[e])}let w;w=0===d||2===d?Math.sqrt(C**2+D**2):Math.sqrt(M**2+F**2);for(let n=f;n<p;n+=g){let i=this.nop[l][n]-1;o(`  - Applied convection boundary condition to node ${i+1} (element ${l+1}, local node ${n+1})`),e[i]+=-s[0]*w*y[n]*u*m;for(let e=f;e<p;e+=g){let o=this.nop[l][e]-1;t[i][o]+=-s[0]*w*y[n]*y[e]*u}}}else if("quadratic"===this.elementOrder)for(let c=0;c<3;c++){let h,f,p,g,b;0===d?(h=n[c],f=0,p=0,g=7,b=3):1===d?(h=0,f=n[c],p=0,g=3,b=1):2===d?(h=n[c],f=1,p=2,g=9,b=3):3===d&&(h=1,f=n[c],p=6,g=9,b=1);let y=a.getBasisFunctions(h,f),v=y.basisFunction,E=y.basisFunctionDerivKsi,C=y.basisFunctionDerivEta,D=0,M=0,F=0,$=0;const w=this.nop[l].length;for(let e=0;e<w;e++){const t=this.nop[l][e]-1;0===d||2===d?(D+=i[t]*E[e],M+=r[t]*E[e]):1!==d&&3!==d||(F+=i[t]*C[e],$+=r[t]*C[e])}let A;A=0===d||2===d?Math.sqrt(D**2+M**2):Math.sqrt(F**2+$**2);for(let n=p;n<g;n+=b){let i=this.nop[l][n]-1;o(`  - Applied convection boundary condition to node ${i+1} (element ${l+1}, local node ${n+1})`),e[i]+=-s[c]*A*v[n]*u*m;for(let e=p;e<g;e+=b){let o=this.nop[l][e]-1;t[i][o]+=-s[c]*A*v[n]*v[e]*u}}}}))}}))}imposeConvectionBoundaryConditionsFront(e,t,n,s,i,r){let a=[],l=[];Object.keys(this.boundaryConditions).forEach((e=>{const t=this.boundaryConditions[e];"convection"===t[0]&&(a[e]=t[1],l[e]=t[2])}));const d=this.nop[e].length,c=Array(d).fill().map((()=>Array(d).fill(0))),u=Array(d).fill(0);for(const m in this.boundaryElements)if("convection"===this.boundaryConditions[m]?.[0]){const h=a[m],f=l[m];o(`Boundary ${m}: Applying convection with heat transfer coefficient h=${h} W/(m²·K) and external temperature T∞=${f} K`);const p=this.boundaryElements[m].find((([t,n])=>t===e));if(p){const a=p[1];if("1D"===this.meshDimension){let t;"linear"===this.elementOrder?t=0===a?0:1:"quadratic"===this.elementOrder&&(t=0===a?0:2),o(`  - Applied convection boundary condition to node ${t+1} (element ${e+1}, local node ${t+1})`),u[t]+=-h*f,c[t][t]+=h}else if("2D"===this.meshDimension)if("linear"===this.elementOrder){let o,l,m,p,g;0===a?(o=s[0],l=0,m=0,p=3,g=2):1===a?(o=0,l=s[0],m=0,p=2,g=1):2===a?(o=s[0],l=1,m=1,p=4,g=2):3===a&&(o=1,l=s[0],m=2,p=4,g=1);const b=r.getBasisFunctions(o,l),y=b.basisFunction,v=b.basisFunctionDerivKsi,E=b.basisFunctionDerivEta;let C,D=0,M=0,F=0,$=0;for(let o=0;o<d;o++){const s=this.nop[e][o]-1;0===a||2===a?(D+=t[s]*v[o],M+=n[s]*v[o]):1!==a&&3!==a||(F+=t[s]*E[o],$+=n[s]*E[o])}C=0===a||2===a?Math.sqrt(D**2+M**2):Math.sqrt(F**2+$**2);for(let e=m;e<p;e+=g){u[e]+=-i[0]*C*y[e]*h*f;for(let t=m;t<p;t+=g)c[e][t]+=-i[0]*C*y[e]*y[t]*h}}else if("quadratic"===this.elementOrder)for(let o=0;o<3;o++){let l,d,m,p,g;0===a?(l=s[o],d=0,m=0,p=7,g=3):1===a?(l=0,d=s[o],m=0,p=3,g=1):2===a?(l=s[o],d=1,m=2,p=9,g=3):3===a&&(l=1,d=s[o],m=6,p=9,g=1);let b=r.getBasisFunctions(l,d),y=b.basisFunction,v=b.basisFunctionDerivKsi,E=b.basisFunctionDerivEta,C=0,D=0,M=0,F=0;const $=this.nop[e].length;for(let o=0;o<$;o++){const s=this.nop[e][o]-1;0===a||2===a?(C+=t[s]*v[o],D+=n[s]*v[o]):1!==a&&3!==a||(M+=t[s]*E[o],F+=n[s]*E[o])}let w;w=0===a||2===a?Math.sqrt(C**2+D**2):Math.sqrt(M**2+F**2);for(let e=m;e<p;e+=g){u[e]+=-i[o]*w*y[e]*h*f;for(let t=m;t<p;t+=g)c[e][t]+=-i[o]*w*y[e]*y[t]*h}}}}return{localJacobianMatrix:c,localResidualVector:u}}}function q(e,t){s("Starting solid heat transfer matrix assembly...");const{nodesXCoordinates:n,nodesYCoordinates:o,nop:i,boundaryElements:r,totalElements:a,meshDimension:l,elementOrder:d}=e,c=X(e),{residualVector:u,jacobianMatrix:m,localToGlobalMap:h,basisFunctions:f,gaussPoints:p,gaussWeights:g,nodesPerElement:b}=c;for(let e=0;e<a;e++){for(let t=0;t<b;t++)h[t]=i[e][t]-1;for(let e=0;e<p.length;e++)if("1D"===l){const t=f.getBasisFunctions(p[e]),o=I({basisFunction:t.basisFunction,basisFunctionDerivKsi:t.basisFunctionDerivKsi,nodesXCoordinates:n,localToGlobalMap:h,nodesPerElement:b}),{detJacobian:s,basisFunctionDerivX:i}=o;for(let t=0;t<b;t++){let n=h[t];for(let o=0;o<b;o++){let r=h[o];m[n][r]+=-g[e]*s*(i[t]*i[o])}}}else if("2D"===l)for(let t=0;t<p.length;t++){const s=f.getBasisFunctions(p[e],p[t]),i=S({basisFunction:s.basisFunction,basisFunctionDerivKsi:s.basisFunctionDerivKsi,basisFunctionDerivEta:s.basisFunctionDerivEta,nodesXCoordinates:n,nodesYCoordinates:o,localToGlobalMap:h,nodesPerElement:b}),{detJacobian:r,basisFunctionDerivX:a,basisFunctionDerivY:l}=i;for(let n=0;n<b;n++){let o=h[n];for(let s=0;s<b;s++){let i=h[s];m[o][i]+=-g[e]*g[t]*r*(a[n]*a[s]+l[n]*l[s])}}}}const y=new j(t,r,i,l,d);return y.imposeConvectionBoundaryConditions(u,m,p,g,n,o,f),y.imposeConstantTempBoundaryConditions(u,m),s("Solid heat transfer matrix assembly completed"),{jacobianMatrix:m,residualVector:u}}function W({elementIndex:e,nop:t,meshData:n,basisFunctions:o,FEAData:s}){const{gaussPoints:i,gaussWeights:r,nodesPerElement:a}=s,{nodesXCoordinates:l,nodesYCoordinates:d,meshDimension:c}=n,u=Array(a).fill().map((()=>Array(a).fill(0))),m=Array(a).fill(0),h=Array(a),f=Array(a);for(let n=0;n<a;n++)h[n]=Math.abs(t[e][n]),f[n]=Math.abs(t[e][n])-1;if("1D"===c)for(let e=0;e<i.length;e++){const{basisFunction:t,basisFunctionDerivKsi:n}=o.getBasisFunctions(i[e]),{detJacobian:s,basisFunctionDerivX:d}=I({basisFunction:t,basisFunctionDerivKsi:n,nodesXCoordinates:l,localToGlobalMap:f,nodesPerElement:a});for(let t=0;t<a;t++)for(let n=0;n<a;n++)u[t][n]-=r[e]*s*(d[t]*d[n])}else if("2D"===c)for(let e=0;e<i.length;e++)for(let t=0;t<i.length;t++){const{basisFunction:n,basisFunctionDerivKsi:s,basisFunctionDerivEta:c}=o.getBasisFunctions(i[e],i[t]),m=h.map((e=>e-1)),{detJacobian:f,basisFunctionDerivX:p,basisFunctionDerivY:g}=S({basisFunction:n,basisFunctionDerivKsi:s,basisFunctionDerivEta:c,nodesXCoordinates:l,nodesYCoordinates:d,localToGlobalMap:m,nodesPerElement:a});for(let n=0;n<a;n++)for(let o=0;o<a;o++)u[n][o]-=r[e]*r[t]*f*(p[n]*p[o]+g[n]*g[o])}return{localJacobianMatrix:u,localResidualVector:m,ngl:h}}class G{constructor(e,t,n,o,s){this.boundaryConditions=e,this.boundaryElements=t,this.nop=n,this.meshDimension=o,this.elementOrder=s}imposeDirichletBoundaryConditions(e,t){"1D"===this.meshDimension?Object.keys(this.boundaryConditions).forEach((n=>{if("constantValue"===this.boundaryConditions[n][0]){const s=this.boundaryConditions[n][1];o(`Boundary ${n}: Applying constant value of ${s} (Dirichlet condition)`),this.boundaryElements[n].forEach((([n,i])=>{if("linear"===this.elementOrder){({0:[0],1:[1]})[i].forEach((i=>{const r=this.nop[n][i]-1;o(`  - Applied constant value to node ${r+1} (element ${n+1}, local node ${i+1})`),e[r]=s;for(let n=0;n<e.length;n++)t[r][n]=0;t[r][r]=1}))}else if("quadratic"===this.elementOrder){({0:[0],1:[2]})[i].forEach((i=>{const r=this.nop[n][i]-1;o(`  - Applied constant value to node ${r+1} (element ${n+1}, local node ${i+1})`),e[r]=s;for(let n=0;n<e.length;n++)t[r][n]=0;t[r][r]=1}))}}))}})):"2D"===this.meshDimension&&Object.keys(this.boundaryConditions).forEach((n=>{if("constantValue"===this.boundaryConditions[n][0]){const s=this.boundaryConditions[n][1];o(`Boundary ${n}: Applying constant value of ${s} (Dirichlet condition)`),this.boundaryElements[n].forEach((([n,i])=>{if("linear"===this.elementOrder){({0:[0,2],1:[0,1],2:[1,3],3:[2,3]})[i].forEach((i=>{const r=this.nop[n][i]-1;o(`  - Applied constant value to node ${r+1} (element ${n+1}, local node ${i+1})`),e[r]=s;for(let n=0;n<e.length;n++)t[r][n]=0;t[r][r]=1}))}else if("quadratic"===this.elementOrder){({0:[0,3,6],1:[0,1,2],2:[2,5,8],3:[6,7,8]})[i].forEach((i=>{const r=this.nop[n][i]-1;o(`  - Applied constant value to node ${r+1} (element ${n+1}, local node ${i+1})`),e[r]=s;for(let n=0;n<e.length;n++)t[r][n]=0;t[r][r]=1}))}}))}}))}imposeConstantValueBoundaryConditionsFront(e,t){"1D"===this.meshDimension?Object.keys(this.boundaryConditions).forEach((n=>{if("constantValue"===this.boundaryConditions[n][0]){const s=this.boundaryConditions[n][1];o(`Boundary ${n}: Applying constant value of ${s} (Dirichlet condition)`),this.boundaryElements[n].forEach((([n,i])=>{if("linear"===this.elementOrder){({0:[0],1:[1]})[i].forEach((i=>{const r=this.nop[n][i]-1;o(`  - Applied constant value to node ${r+1} (element ${n+1}, local node ${i+1})`),e[r]=1,t[r]=s}))}else if("quadratic"===this.elementOrder){({0:[0],2:[2]})[i].forEach((i=>{const r=this.nop[n][i]-1;o(`  - Applied constant value to node ${r+1} (element ${n+1}, local node ${i+1})`),e[r]=1,t[r]=s}))}}))}})):"2D"===this.meshDimension&&Object.keys(this.boundaryConditions).forEach((n=>{if("constantValue"===this.boundaryConditions[n][0]){const s=this.boundaryConditions[n][1];o(`Boundary ${n}: Applying constant value of ${s} (Dirichlet condition)`),this.boundaryElements[n].forEach((([n,i])=>{if("linear"===this.elementOrder){({0:[0,2],1:[0,1],2:[1,3],3:[2,3]})[i].forEach((i=>{const r=this.nop[n][i]-1;o(`  - Applied constant value to node ${r+1} (element ${n+1}, local node ${i+1})`),e[r]=1,t[r]=s}))}else if("quadratic"===this.elementOrder){({0:[0,3,6],1:[0,1,2],2:[2,5,8],3:[6,7,8]})[i].forEach((i=>{const r=this.nop[n][i]-1;o(`  - Applied constant value to node ${r+1} (element ${n+1}, local node ${i+1})`),e[r]=1,t[r]=s}))}}))}}))}}function K(e,t,n,i){s("Starting front propagation matrix assembly...");let r=1-i+.01;o(`eikonalViscousTerm: ${r}`),o(`eikonalActivationFlag: ${i}`);const{nodesXCoordinates:a,nodesYCoordinates:l,nop:d,boundaryElements:c,totalElements:u,meshDimension:m,elementOrder:h}=e,f=X(e),{residualVector:p,jacobianMatrix:g,localToGlobalMap:b,basisFunctions:y,gaussPoints:v,gaussWeights:E,nodesPerElement:C}=f;for(let e=0;e<u;e++){for(let t=0;t<C;t++)b[t]=d[e][t]-1;for(let e=0;e<v.length;e++)if("1D"===m){errorLog("1D front propagation is not yet supported");let t=y.getBasisFunctions(v[e]);const o=I({basisFunction:t.basisFunction,basisFunctionDerivKsi:t.basisFunctionDerivKsi,nodesXCoordinates:a,localToGlobalMap:b,nodesPerElement:C}),{detJacobian:s,basisFunctionDerivX:i}=o;t.basisFunction;let r=0;for(let e=0;e<C;e++)r+=n[b[e]]*i[e];for(let e=0;e<C;e++){b[e];for(let e=0;e<C;e++)b[e]}}else if("2D"===m)for(let t=0;t<v.length;t++){let o=y.getBasisFunctions(v[e],v[t]);const s=S({basisFunction:o.basisFunction,basisFunctionDerivKsi:o.basisFunctionDerivKsi,basisFunctionDerivEta:o.basisFunctionDerivEta,nodesXCoordinates:a,nodesYCoordinates:l,localToGlobalMap:b,nodesPerElement:C}),{detJacobian:d,basisFunctionDerivX:c,basisFunctionDerivY:u}=s,m=o.basisFunction;let h=0,f=0;for(let e=0;e<C;e++)h+=n[b[e]]*c[e],f+=n[b[e]]*u[e];for(let n=0;n<C;n++){let o=b[n];p[o]+=r*E[e]*E[t]*d*c[n]*h+r*E[e]*E[t]*d*u[n]*f,0!==i&&(p[o]+=i*(E[e]*E[t]*d*m[n]*Math.sqrt(h**2+f**2)-E[e]*E[t]*d*m[n]));for(let s=0;s<C;s++){let a=b[s];g[o][a]+=-r*E[e]*E[t]*d*(c[n]*c[s]+u[n]*u[s]),0!==i&&(g[o][a]+=i*(-d*h*m[n]*E[e]*E[t]/Math.sqrt(h**2+f**2+1e-8))*c[s]-i*(d*f*m[n]*E[e]*E[t]/Math.sqrt(h**2+f**2+1e-8))*u[s])}}}}return new G(t,c,d,m,h).imposeDirichletBoundaryConditions(p,g),s("Front propagation matrix assembly completed"),{jacobianMatrix:g,residualVector:p}}function L({elementIndex:e,nop:t,meshData:n,basisFunctions:o,FEAData:s,solutionVector:i,eikonalActivationFlag:r}){const{gaussPoints:a,gaussWeights:l,nodesPerElement:d}=s,{nodesXCoordinates:c,nodesYCoordinates:u,meshDimension:m}=n;let h=1-r+.01;const f=Array(d).fill().map((()=>Array(d).fill(0))),p=Array(d).fill(0),g=Array(d),b=Array(d);for(let n=0;n<d;n++)g[n]=Math.abs(t[e][n]),b[n]=Math.abs(t[e][n])-1;for(let e=0;e<a.length;e++)if("1D"===m){errorLog("1D front propagation is not yet supported");let t=o.getBasisFunctions(a[e]);const n=I({basisFunction:t.basisFunction,basisFunctionDerivKsi:t.basisFunctionDerivKsi,nodesXCoordinates:c,localToGlobalMap:b,nodesPerElement:d}),{detJacobian:s,basisFunctionDerivX:r}=n;t.basisFunction;let l=0;for(let e=0;e<d;e++)l+=i[b[e]]*r[e];for(let e=0;e<d;e++){b[e];for(let e=0;e<d;e++)b[e]}}else if("2D"===m)for(let t=0;t<a.length;t++){const{basisFunction:n,basisFunctionDerivKsi:s,basisFunctionDerivEta:m}=o.getBasisFunctions(a[e],a[t]),{detJacobian:g,basisFunctionDerivX:y,basisFunctionDerivY:v}=S({basisFunction:n,basisFunctionDerivKsi:s,basisFunctionDerivEta:m,nodesXCoordinates:c,nodesYCoordinates:u,localToGlobalMap:b,nodesPerElement:d});let E=0,C=0;for(let e=0;e<d;e++)E+=i[b[e]]*y[e],C+=i[b[e]]*v[e];for(let o=0;o<d;o++){b[o],p[o]+=h*l[e]*l[t]*g*y[o]*E+h*l[e]*l[t]*g*v[o]*C,0!==r&&(p[o]+=r*(l[e]*l[t]*g*n[o]*Math.sqrt(E**2+C**2)-l[e]*l[t]*g*n[o]));for(let s=0;s<d;s++)f[o][s]-=h*l[e]*l[t]*g*(y[o]*y[s]+v[o]*v[s]),0!==r&&(f[o][s]+=r*(-g*E*n[o]*l[e]*l[t]/Math.sqrt(E**2+C**2+1e-8))*y[s]-r*(g*C*n[o]*l[e]*l[t]/Math.sqrt(E**2+C**2+1e-8))*v[s])}}return{localJacobianMatrix:f,localResidualVector:p,ngl:g}}const J={},H={},U={currentElementIndex:0},z={};let _;function Q(e,t,n,r={}){const a=X(t),l=t.nodesXCoordinates.length,d=t.totalElements;!function(e,t){J.nodalNumbering=Array(t).fill().map((()=>Array(e).fill(0))),J.nodeConstraintCode=Array(e).fill(0),J.boundaryValues=Array(e).fill(0),J.globalResidualVector=Array(e).fill(0),J.solutionVector=Array(e).fill(0),J.topologyData=Array(t).fill(0),J.lateralData=Array(t).fill(0),H.writeFlag=0,H.totalNodes=e,H.transformationFlag=0,H.nodesPerElement=Array(t).fill(0),H.determinant=1;const n=Math.max(e,2e3);H.globalSolutionVector=Array(n).fill(0),H.frontDataIndex=0,U.localJacobianMatrix=Array(e).fill().map((()=>Array(e).fill(0))),U.currentElementIndex=0;const o=function(e,t){const n=Math.max(Math.ceil(Math.sqrt(t))*e,2*e);return n*t}(e,t);z.frontValues=Array(o).fill(0),z.columnHeaders=Array(n).fill(0),z.pivotRow=Array(n).fill(0),z.pivotData=Array(o).fill(0)}(a.nodesPerElement,d),s("Solving system using frontal..."),console.time("systemSolving"),_=new x({meshDimension:t.meshDimension,elementOrder:t.elementOrder});for(let e=0;e<t.totalElements;e++)for(let n=0;n<a.nodesPerElement;n++)J.nodalNumbering[e][n]=t.nop[e][n];for(let e=0;e<t.nodesXCoordinates.length;e++)J.nodeConstraintCode[e]=0,J.boundaryValues[e]=0;let c;e===W?(c=new j(n,t.boundaryElements,t.nop,t.meshDimension,t.elementOrder),c.imposeConstantTempBoundaryConditionsFront(J.nodeConstraintCode,J.boundaryValues)):e===L&&(c=new G(n,t.boundaryElements,t.nop,t.meshDimension,t.elementOrder),c.imposeConstantValueBoundaryConditionsFront(J.nodeConstraintCode,J.boundaryValues));for(let e=0;e<t.nodesXCoordinates.length;e++)J.globalResidualVector[e]=0;H.totalNodes=t.nodesXCoordinates.length,H.writeFlag=0,H.transformationFlag=1,H.determinant=1;for(let e=0;e<t.totalElements;e++)H.nodesPerElement[e]=a.nodesPerElement;H.currentSolutionVector=r.solutionVector,H.eikonalActivationFlag=r.eikonalActivationFlag,function(e,t,n,s){const r=e.totalElements,a=e.nodesXCoordinates.length,l=Math.max(a,H.globalSolutionVector.length);let d,c=Array(t.nodesPerElement).fill(0),u=Array(t.nodesPerElement).fill(0),m=Array(l).fill(0),h=Array(l).fill(0),f=Array(l).fill(0),p=Array(l).fill(0),g=Array(l).fill(0),b=Array(l).fill().map((()=>Array(l).fill(0))),y=Array(a).fill(0),v=Array(a).fill(0),E=Array(a).fill(0),C=1;H.writeFlag++;let D=1,M=1;U.currentElementIndex=0;for(let e=0;e<H.totalNodes;e++)y[e]=0,v[e]=0;if(0!==H.transformationFlag){for(let e=0;e<H.totalNodes;e++)E[e]=0;for(let e=0;e<r;e++){let t=r-e-1;for(let e=0;e<H.nodesPerElement[t];e++){let n=J.nodalNumbering[t][e];0===E[n-1]&&(E[n-1]=1,J.nodalNumbering[t][e]=-J.nodalNumbering[t][e])}}}H.transformationFlag=0;let F=0,$=0;for(let e=0;e<l;e++)for(let t=0;t<l;t++)b[t][e]=0;for(;;){let a=!1,E=0,w=0;if(U.currentElementIndex<r&&(U.currentElementIndex++,a=Z(e,t,n,s)),a){const e=U.currentElementIndex;E=H.nodesPerElement[e-1],w=H.nodesPerElement[e-1];for(let t=0;t<w;t++){let n,o,s=J.nodalNumbering[e-1][t];if(0===F)F++,c[t]=F,z.columnHeaders[F-1]=s;else{for(n=0;n<F&&Math.abs(s)!==Math.abs(z.columnHeaders[n]);n++);n===F?(F++,c[t]=F,z.columnHeaders[F-1]=s):(c[t]=n+1,z.columnHeaders[n]=s)}if(0===$)$++,u[t]=$,m[$-1]=s;else{for(o=0;o<$&&Math.abs(s)!==Math.abs(m[o]);o++);o===$?($++,u[t]=$,m[$-1]=s):(u[t]=o+1,m[o]=s)}}if($>l||F>l)return void i("Error: systemSize not large enough");for(let e=0;e<w;e++){let t=c[e];for(let n=0;n<E;n++){b[u[n]-1][t-1]+=U.localJacobianMatrix[n][e]}}}let A=0;for(let e=0;e<F;e++)z.columnHeaders[e]<0&&(f[A]=e+1,A++);let x=0,N=0;for(let e=0;e<$;e++){let t=m[e];if(t<0){h[N]=e+1,N++;let n=Math.abs(t);1===J.nodeConstraintCode[n-1]&&(p[x]=e+1,x++,J.nodeConstraintCode[n-1]=2,J.globalResidualVector[n-1]=J.boundaryValues[n-1])}}if(x>0)for(let e=0;e<x;e++){let t=p[e]-1,n=Math.abs(m[t]);for(let e=0;e<F;e++){b[t][e]=0,Math.abs(z.columnHeaders[e])===n&&(b[t][e]=1)}}if(A>M||U.currentElementIndex<r){if(0===A)return void i("Error: no more rows fully summed");let e=h[0],t=f[0],n=b[e-1][t-1];if(Math.abs(n)<1e-4){n=0;for(let o=0;o<A;o++){let s=f[o];for(let o=0;o<N;o++){let i=h[o],r=b[i-1][s-1];Math.abs(r)>Math.abs(n)&&(n=r,t=s,e=i)}}}let s=Math.abs(m[e-1]);d=Math.abs(z.columnHeaders[t-1]);let a=s+d+y[s-1]+v[d-1];H.determinant=H.determinant*n*(-1)**a/Math.abs(n);for(let e=0;e<H.totalNodes;e++)e>=s&&y[e]--,e>=d&&v[e]--;if(Math.abs(n)<1e-10&&i(`Matrix singular or ill-conditioned, currentElementIndex=${U.currentElementIndex}, pivotGlobalRowIndex=${s}, pivotColumnGlobalIndex=${d}, pivotValue=${n}`),0===n)return;for(let t=0;t<F;t++)z.pivotRow[t]=b[e-1][t]/n;let l=J.globalResidualVector[s-1]/n;if(J.globalResidualVector[s-1]=l,g[e-1]=n,e>1)for(let n=0;n<e-1;n++){let e=Math.abs(m[n]),o=b[n][t-1];if(g[n]=o,t>1&&0!==o)for(let e=0;e<t-1;e++)b[n][e]-=o*z.pivotRow[e];if(t<F)for(let e=t;e<F;e++)b[n][e-1]=b[n][e]-o*z.pivotRow[e];J.globalResidualVector[e-1]-=o*l}if(e<$)for(let n=e;n<$;n++){let e=Math.abs(m[n]),o=b[n][t-1];if(g[n]=o,t>1)for(let e=0;e<t-1;e++)b[n-1][e]=b[n][e]-o*z.pivotRow[e];if(t<F)for(let e=t;e<F;e++)b[n-1][e-1]=b[n][e]-o*z.pivotRow[e];J.globalResidualVector[e-1]-=o*l}for(let e=0;e<$;e++)z.pivotData[D+e-1]=g[e];D+=$;for(let e=0;e<$;e++)z.pivotData[D+e-1]=m[e];D+=$,z.pivotData[D-1]=e,D++;for(let e=0;e<F;e++)z.frontValues[C-1+e]=z.pivotRow[e];C+=F;for(let e=0;e<F;e++)z.frontValues[C-1+e]=z.columnHeaders[e];C+=F,z.frontValues[C-1]=s,z.frontValues[C]=F,z.frontValues[C+1]=t,z.frontValues[C+2]=n,C+=4;for(let e=0;e<$;e++)b[e][F-1]=0;for(let e=0;e<F;e++)b[$-1][e]=0;if(F--,t<F+1)for(let e=t-1;e<F;e++)z.columnHeaders[e]=z.columnHeaders[e+1];if($--,e<$+1)for(let t=e-1;t<$;t++)m[t]=m[t+1];if($>1||U.currentElementIndex<r)continue;if(d=Math.abs(z.columnHeaders[0]),e=1,n=b[0][0],s=Math.abs(m[0]),t=1,a=s+d+y[s-1]+v[d-1],H.determinant=H.determinant*n*(-1)**a/Math.abs(n),z.pivotRow[0]=1,Math.abs(n)<1e-10&&i(`Matrix singular or ill-conditioned, currentElementIndex=${U.currentElementIndex}, pivotGlobalRowIndex=${s}, pivotColumnGlobalIndex=${d}, pivotValue=${n}`),0===n)return;J.globalResidualVector[s-1]=J.globalResidualVector[s-1]/n,z.frontValues[C-1]=z.pivotRow[0],C++,z.frontValues[C-1]=z.columnHeaders[0],C++,z.frontValues[C-1]=s,z.frontValues[C]=F,z.frontValues[C+1]=t,z.frontValues[C+2]=n,C+=4,z.pivotData[D-1]=g[0],D++,z.pivotData[D-1]=m[0],D++,z.pivotData[D-1]=e,D++,H.frontDataIndex=C,1===H.writeFlag&&o(`total ecs transfer in matrix reduction=${C}`),ee(C);break}}}(t,a,c,e);for(let e=0;e<t.nodesXCoordinates.length;e++)J.solutionVector[e]=H.globalSolutionVector[e];const{nodesXCoordinates:u,nodesYCoordinates:m}=t;for(let e=0;e<t.nodesXCoordinates.length;e++)"1D"===t.meshDimension?o(`${u[e].toExponential(5)}  ${J.solutionVector[e].toExponential(5)}`):o(`${u[e].toExponential(5)}  ${m[e].toExponential(5)}  ${J.solutionVector[e].toExponential(5)}`);console.timeEnd("systemSolving"),s("System solved successfully");const{nodesXCoordinates:h,nodesYCoordinates:f}=t;return{solutionVector:J.solutionVector.slice(0,l),nodesCoordinates:{nodesXCoordinates:h,nodesYCoordinates:f}}}function Z(e,t,n,o){const s=U.currentElementIndex-1;if(s<0||s>=e.totalElements)return i(`Skipping out-of-range elementIndex=${s} (totalElements=${e.totalElements})`),!1;const{localJacobianMatrix:r,localResidualVector:a,ngl:l}=o({elementIndex:s,nop:J.nodalNumbering,meshData:e,basisFunctions:_,FEAData:t,solutionVector:H.currentSolutionVector,eikonalActivationFlag:H.eikonalActivationFlag});let d=Array(t.nodesPerElement).fill().map((()=>Array(t.nodesPerElement).fill(0))),c=Array(t.nodesPerElement).fill(0);if(o===W){let o=!1;for(const t in e.boundaryElements)if("convection"===n.boundaryConditions[t]?.[0]&&e.boundaryElements[t].some((([e,t])=>e===s))){o=!0;break}if(o){const{gaussPoints:o,gaussWeights:i}=t,r=n.imposeConvectionBoundaryConditionsFront(s,e.nodesXCoordinates,e.nodesYCoordinates,o,i,_);d=r.localJacobianMatrix,c=r.localResidualVector}}for(let e=0;e<t.nodesPerElement;e++)for(let n=0;n<t.nodesPerElement;n++)U.localJacobianMatrix[e][n]=r[e][n]+d[e][n];for(let e=0;e<t.nodesPerElement;e++){const t=l[e]-1;J.globalResidualVector[t]+=a[e]+c[e]}return!0}function ee(e){for(let e=0;e<H.totalNodes;e++)H.globalSolutionVector[e]=J.boundaryValues[e];for(let t=1;t<=H.totalNodes;t++){e-=4;let n=z.frontValues[e-1],o=z.frontValues[e],s=z.frontValues[e+1];if(z.frontValues[e+2],1===t)e--,z.columnHeaders[0]=z.frontValues[e-1],e--,z.pivotRow[0]=z.frontValues[e-1];else{e-=o;for(let t=0;t<o;t++)z.columnHeaders[t]=z.frontValues[e-1+t];e-=o;for(let t=0;t<o;t++)z.pivotRow[t]=z.frontValues[e-1+t]}let i=Math.abs(z.columnHeaders[s-1]);if(J.nodeConstraintCode[i-1]>0)continue;let r=0;z.pivotRow[s-1]=0;for(let e=0;e<o;e++)r-=z.pivotRow[e]*H.globalSolutionVector[Math.abs(z.columnHeaders[e])-1];H.globalSolutionVector[i-1]=r+J.globalResidualVector[n-1],J.nodeConstraintCode[i-1]=1}1===H.writeFlag&&o(`value of frontDataCounter after backsubstitution=${e}`)}function te(t,n={}){let o=0,r=!1,a=0,l=[],d=[],c=[],u=[];const{maxIterations:m=100,tolerance:h=1e-4}=n;let f=n.meshData.nodesXCoordinates.length;for(let e=0;e<f;e++)l[e]=0,d[e]=0;for(n.initialSolution&&n.initialSolution.length===f&&(d=[...n.initialSolution]);a<m&&!r;){for(let e=0;e<d.length;e++)d[e]=Number(d[e])+Number(l[e]);if("frontal"===n.solverMethod){l=Q(L,n.meshData,n.boundaryConditions,{solutionVector:d,eikonalActivationFlag:n.eikonalActivationFlag}).solutionVector}else{({jacobianMatrix:c,residualVector:u}=t(n.meshData,n.boundaryConditions,d,n.eikonalActivationFlag));l=w(n.solverMethod,c,u).solutionVector}if(o=e(l),s(`Newton-Raphson iteration ${a+1}: Error norm = ${o.toExponential(4)}`),o<=h)r=!0;else if(o>100){i(`Solution not converged. Error norm: ${o}`);break}a++}return{solutionVector:d,converged:r,iterations:a,jacobianMatrix:c,residualVector:u}}class ne{constructor(e,t,n,o,s,i,r,a){this.boundaryConditions=e,this.boundaryElements=t,this.nop=n,this.meshDimension=o,this.elementOrder=s,this.totalNodesVelocity=i,this.totalNodesPressure=r,this.q2ToPressureMap=a}imposeDirichletBoundaryConditions(e,t){const n=e.length;let s=!1;if("2D"===this.meshDimension&&(Object.keys(this.boundaryConditions).forEach((i=>{const r=this.boundaryConditions[i][0];if("stressFree"===r)s=!0,o(`Boundary ${i}: Applying stress-free condition (natural BC)`);else if("constantVelocity"===r){const s=this.boundaryConditions[i][1],r=this.boundaryConditions[i][2];o(`Boundary ${i}: Applying constant velocity condition (u=${s}, v=${r})`),this.boundaryElements[i].forEach((([i,a])=>{if("quadratic"===this.elementOrder){({0:[0,3,6],1:[0,1,2],2:[2,5,8],3:[6,7,8]})[a].forEach((a=>{const l=this.nop[i][a]-1,d=l,c=this.totalNodesVelocity+l;o(`  - Applied velocity Dirichlet to node ${l+1} (element ${i+1}, local node ${a+1})`),e[d]=s;for(let e=0;e<n;e++)t[d][e]=0;t[d][d]=1,e[c]=r;for(let e=0;e<n;e++)t[c][e]=0;t[c][c]=1}))}else if("linear"===this.elementOrder){({0:[0,2],1:[0,1],2:[1,3],3:[2,3]})[a].forEach((a=>{const l=this.nop[i][a]-1,d=l,c=this.totalNodesVelocity+l;o(`  - Applied velocity Dirichlet to node ${l+1} (element ${i+1}, local node ${a+1})`),e[d]=s;for(let e=0;e<n;e++)t[d][e]=0;t[d][d]=1,e[c]=r;for(let e=0;e<n;e++)t[c][e]=0;t[c][c]=1}))}}))}})),!s)){const s=2*this.totalNodesVelocity;for(let e=0;e<n;e++)t[s][e]=0;t[s][s]=1,e[s]=0,o("Pinned pressure at first pressure node (p = 0) to remove null space")}}}class oe{constructor(){var e;this.solverConfig=null,this.meshConfig={},this.boundaryConditions={},this.solverMethod="lusolve",this.coefficientFunctions=null,e="FEAScript is provided “as is” without any warranty. The authors are not responsible for any damages or losses that may result from using the software. See the license for more details: https://github.com/FEAScript/FEAScript-core/blob/main/LICENSE",console.log("%c[WARN] "+e,"color: #FF9800; font-weight: bold;"),s("FEAScriptModel instance created")}setModelConfig(e,t={}){this.solverConfig=e,void 0!==t?.coefficientFunctions&&(this.coefficientFunctions=t.coefficientFunctions,o("coefficientFunctions set")),void 0!==t?.maxIterations&&(this.maxIterations=t.maxIterations,o(`maxIterations set to ${this.maxIterations}`)),void 0!==t?.tolerance&&(this.tolerance=t.tolerance,o(`tolerance set to ${this.tolerance}`)),o(`solverConfig set to ${e}`)}setMeshConfig(e){this.meshConfig=e,o(`meshConfig set with dimensions: ${e.meshDimension}`)}addBoundaryCondition(e,t){this.boundaryConditions[e]=t,o(`boundaryConditions added for boundary: ${e}, type: ${t[0]}`)}setSolverMethod(e){this.solverMethod=e,o(`solverMethod set to: ${e}`)}solve(e={}){this.solverConfig&&this.meshConfig&&this.boundaryConditions||i("solverConfig, meshConfig and boundaryConditions must be set before solving");let t=[],n=[],r=[],a=[];s("Preparing mesh...");const l=V(this.meshConfig);s("Mesh preparation completed");const d={nodesXCoordinates:l.nodesXCoordinates,nodesYCoordinates:l.nodesYCoordinates};if(s("Beginning solving process..."),console.time("totalSolvingTime"),s(`Using solver ${this.solverConfig}`),"heatConductionScript"===this.solverConfig)if("frontal"===this.solverMethod){r=Q(W,l,this.boundaryConditions).solutionVector}else{({jacobianMatrix:t,residualVector:n}=q(l,this.boundaryConditions));r=w(this.solverMethod,t,n,{maxIterations:e.maxIterations??this.maxIterations,tolerance:e.tolerance??this.tolerance}).solutionVector}else if("frontPropagationScript"===this.solverConfig){let o=0;const s=5,i={meshData:l,boundaryConditions:this.boundaryConditions,eikonalActivationFlag:o,solverMethod:this.solverMethod,initialSolution:a,maxIterations:e.maxIterations??this.maxIterations,tolerance:e.tolerance??this.tolerance};for(;o<=1;){i.eikonalActivationFlag=o,r.length>0&&(i.initialSolution=[...r]);const e=te(K,i);t=e.jacobianMatrix,n=e.residualVector,r=e.solutionVector,o+=1/s}}else if("generalFormPDEScript"===this.solverConfig)if("frontal"===this.solverMethod)i("Frontal solver is not yet supported for generalFormPDEScript. Please use 'lusolve' or 'jacobi'.");else{({jacobianMatrix:t,residualVector:n}=function(e,t,n){s("Starting general form PDE matrix assembly...");const{nodesXCoordinates:o,nodesYCoordinates:r,nop:a,boundaryElements:l,totalElements:d,meshDimension:c,elementOrder:u}=e,{A:m,B:h,C:f,D:p}=n,g=X(e),{residualVector:b,jacobianMatrix:y,localToGlobalMap:v,basisFunctions:E,gaussPoints:C,gaussWeights:D,nodesPerElement:M}=g;if("1D"===c)for(let e=0;e<d;e++){for(let t=0;t<M;t++)v[t]=Math.abs(a[e][t])-1;for(let e=0;e<C.length;e++){const{basisFunction:t,basisFunctionDerivKsi:n}=E.getBasisFunctions(C[e]),{detJacobian:s,basisFunctionDerivX:i}=I({basisFunction:t,basisFunctionDerivKsi:n,nodesXCoordinates:o,localToGlobalMap:v,nodesPerElement:M});let r=0;for(let e=0;e<M;e++)r+=o[v[e]]*t[e];const a=m(r),l=h(r),d=f(r),c=p(r);for(let n=0;n<M;n++){const o=v[n];b[o]-=D[e]*s*c*t[n];for(let r=0;r<M;r++){const c=v[r];y[o][c]+=D[e]*s*a*i[n]*i[r],y[o][c]-=D[e]*s*l*i[r]*t[n],y[o][c]-=D[e]*s*d*t[n]*t[r]}}}}else"2D"===c&&i("2D general form PDE is not yet supported in assembleGeneralFormPDEMat.");return new G(t,l,a,c,u).imposeDirichletBoundaryConditions(b,y),s("General form PDE matrix assembly completed"),{jacobianMatrix:y,residualVector:b}}(l,this.boundaryConditions,this.coefficientFunctions));r=w(this.solverMethod,t,n,{maxIterations:e.maxIterations??this.maxIterations,tolerance:e.tolerance??this.tolerance}).solutionVector}else if("creepingFlowScript"===this.solverConfig){const a=function(e,t){s("Starting creeping flow matrix assembly...");const{nodesXCoordinates:n,nodesYCoordinates:r,nop:a,boundaryElements:l,totalElements:d,totalNodes:c,meshDimension:u,elementOrder:m}=e;"2D"!==u&&i("Creeping flow solver requires a 2D mesh"),"quadratic"!==m&&i("Creeping flow solver requires quadratic elements for Taylor-Hood (Q2-Q1) formulation");const h=c,f=[0,2,6,8],p=new Map,g=[];let b=0;for(let e=0;e<d;e++)for(let t=0;t<f.length;t++){const n=a[e][f[t]]-1;p.has(n)||(p.set(n,b),g.push(n),b++)}const y=b,v=2*h+y;o(`Creeping flow DOFs: ${h} velocity nodes (Q2), ${y} pressure nodes (Q1), ${v} total DOFs`);let E=[],C=[];for(let e=0;e<v;e++){E[e]=0,C.push([]);for(let t=0;t<v;t++)C[e][t]=0}const D=new x({meshDimension:"2D",elementOrder:"quadratic"}),M=new x({meshDimension:"2D",elementOrder:"linear"});let F=new O({meshDimension:"2D",elementOrder:"quadratic"}).getGaussPointsAndWeights(),$=F.gaussPoints,w=F.gaussWeights;for(let e=0;e<d;e++){let t=[];for(let n=0;n<9;n++)t[n]=a[e][n]-1;let o=[];for(let t=0;t<4;t++){const n=a[e][f[t]]-1;o[t]=p.get(n)}for(let e=0;e<$.length;e++)for(let s=0;s<$.length;s++){const i=D.getBasisFunctions($[e],$[s]),a=M.getBasisFunctions($[e],$[s]),l=S({basisFunction:i.basisFunction,basisFunctionDerivKsi:i.basisFunctionDerivKsi,basisFunctionDerivEta:i.basisFunctionDerivEta,nodesXCoordinates:n,nodesYCoordinates:r,localToGlobalMap:t,nodesPerElement:9}),{detJacobian:d,basisFunctionDerivX:c,basisFunctionDerivY:u}=l,m=w[e]*w[s]*d;for(let e=0;e<9;e++){let n=t[e],s=n,i=h+n;for(let n=0;n<9;n++){let o=t[n],r=o,a=h+o,l=1*-m*(c[e]*c[n]+u[e]*u[n]);C[s][r]+=l,C[i][a]+=l}for(let t=0;t<4;t++){let n=2*h+o[t],r=m*a.basisFunction[t]*c[e],l=m*a.basisFunction[t]*u[e];C[s][n]+=r,C[i][n]+=l,C[n][s]+=-r,C[n][i]+=-l}}}}return new ne(t,l,a,u,m,h,y,p).imposeDirichletBoundaryConditions(E,C),s("Creeping flow matrix assembly completed"),{jacobianMatrix:C,residualVector:E,totalNodesVelocity:h,totalNodesPressure:y,pressureNodeIndices:g}}(l,this.boundaryConditions);t=a.jacobianMatrix,n=a.residualVector;r=w(this.solverMethod,t,n,{maxIterations:e.maxIterations??this.maxIterations,tolerance:e.tolerance??this.tolerance}).solutionVector,this._creepingFlowMetadata={totalNodesVelocity:a.totalNodesVelocity,totalNodesPressure:a.totalNodesPressure,pressureNodeIndices:a.pressureNodeIndices}}return console.timeEnd("totalSolvingTime"),s("Solving process completed"),{solutionVector:r,nodesCoordinates:d}}async solveAsync(e,t={}){this.solverConfig&&this.meshConfig&&this.boundaryConditions||i("Solver config, mesh config, and boundary conditions must be set before solving.");let n=[],o=[],r=[];s("Preparing mesh...");const a=V(this.meshConfig);s("Mesh preparation completed");const l={nodesXCoordinates:a.nodesXCoordinates,nodesYCoordinates:a.nodesYCoordinates};if(s("Beginning solving process..."),console.time("totalSolvingTime"),s(`Using solver: ${this.solverConfig}`),"heatConductionScript"===this.solverConfig&&(({jacobianMatrix:n,residualVector:o}=q(a,this.boundaryConditions)),"jacobi-gpu"===this.solverMethod)){const{solutionVector:s}=await A("jacobi-gpu",n,o,{computeEngine:e,maxIterations:t.maxIterations??this.maxIterations,tolerance:t.tolerance??this.tolerance});r=s}return console.timeEnd("totalSolvingTime"),s("Solving process completed"),{solutionVector:r,nodesCoordinates:l}}}const se=async e=>{let t={nodesXCoordinates:[],nodesYCoordinates:[],nodalNumbering:{quadElements:[],triangleElements:[]},boundaryElements:[],boundaryConditions:[],boundaryNodePairs:{},gmshV:0,ascii:!1,fltBytes:"8",totalNodesX:0,totalNodesY:0,physicalPropMap:[],elementTypes:{}};const n={curves:{}};let s=(await e.text()).split("\n").map((e=>e.trim())).filter((e=>""!==e)),i="",r=0,a=0,l=0,d=0,c={numNodes:0},u=[],m=0,h=0,f=0,p=0,g={numElements:0},b=0,y={},v=null,E=null,C=0,D=0,M=0;for(;r<s.length;){const e=s[r];if("$MeshFormat"===e){i="meshFormat",r++;continue}if("$EndMeshFormat"===e){i="",r++;continue}if("$PhysicalNames"===e){i="physicalNames",r++;continue}if("$EndPhysicalNames"===e){i="",r++;continue}if("$Entities"===e){i="entities",E="counts",r++;continue}if("$EndEntities"===e){i="",v=null,E=null,r++;continue}if("$Nodes"===e){i="nodes",r++;continue}if("$EndNodes"===e){i="",r++;continue}if("$Elements"===e){i="elements",r++;continue}if("$EndElements"===e){i="",r++;continue}const o=e.split(/\s+/);if("meshFormat"===i)t.gmshV=parseFloat(o[0]),t.ascii="0"===o[1],t.fltBytes=o[2];else if("physicalNames"===i){const e=parseInt(o[0],10),n=parseInt(o[1],10);let s=o.slice(2).join(" ").replace(/^"|"$/g,"");t.physicalPropMap.push({tag:n,dimension:e,name:s})}else if("entities"===i)if("counts"===E)v={points:parseInt(o[0],10),curves:parseInt(o[1],10),surfaces:parseInt(o[2],10),volumes:parseInt(o[3],10)},E="points";else if("points"===E)C++,C===v.points&&(E="curves");else if("curves"===E){const e=parseInt(o[0],10),t=parseInt(o[7],10),s=o.slice(8,8+t).map((e=>parseInt(e,10)));n.curves[e]=s,D++,D===v.curves&&(E="surfaces")}else"surfaces"===E&&(M++,M===v.surfaces&&(E="volumes"));else if("nodes"===i){if(0===a)a=parseInt(o[0],10),l=parseInt(o[1],10),t.nodesXCoordinates=new Array(l).fill(0),t.nodesYCoordinates=new Array(l).fill(0);else if(d<a&&0===c.numNodes)c={dim:parseInt(o[0],10),tag:parseInt(o[1],10),parametric:parseInt(o[2],10),numNodes:parseInt(o[3],10)},u=[],m=0,h=0;else if(m<c.numNodes){for(let e of o)if(u.push(parseInt(e,10)),m++,m===c.numNodes)break}else if(h<c.numNodes){const e=u[h]-1;t.nodesXCoordinates[e]=parseFloat(o[0]),t.nodesYCoordinates[e]=parseFloat(o[1]),t.totalNodesX++,t.totalNodesY++,h++,h===c.numNodes&&(d++,c={numNodes:0})}}else if("elements"===i)if(0===f)f=parseInt(o[0],10),parseInt(o[1],10);else if(p<f&&0===g.numElements)g={dim:parseInt(o[0],10),tag:parseInt(o[1],10),elementType:parseInt(o[2],10),numElements:parseInt(o[3],10)},t.elementTypes[g.elementType]=(t.elementTypes[g.elementType]||0)+g.numElements,b=0;else if(b<g.numElements){const e=o.slice(1).map((e=>parseInt(e,10)));let s=g.tag;if(1===g.dim){const e=n.curves[g.tag];e&&e.length>0&&(s=e[0])}1===g.elementType||8===g.elementType?(y[s]||(y[s]=[]),y[s].push(e),t.boundaryNodePairs[s]||(t.boundaryNodePairs[s]=[]),t.boundaryNodePairs[s].push(e)):2===g.elementType?t.nodalNumbering.triangleElements.push(e):3!==g.elementType&&10!==g.elementType||t.nodalNumbering.quadElements.push(e),b++,b===g.numElements&&(p++,g={numElements:0})}r++}return t.physicalPropMap.forEach((e=>{if(1===e.dimension){const n=y[e.tag]||[];n.length>0&&t.boundaryConditions.push({name:e.name,tag:e.tag,nodes:n})}})),o(`Parsed boundary node pairs by physical tag: ${JSON.stringify(t.boundaryNodePairs)}`),t};function ie(e,t,n,o){console.time("plottingTime");const{nodesXCoordinates:s,nodesYCoordinates:i}=t.nodesCoordinates,r=t.solutionVector,a=e.solverConfig,l=e.meshConfig.meshDimension;if(V(e.meshConfig),"1D"===l&&"line"===n){let e;e=r.length>0&&Array.isArray(r[0])?r.map((e=>e[0])):r,Array.from(s);let t={x:s,y:e,mode:"lines",type:"scatter",line:{color:"rgb(219, 64, 82)",width:2},name:"Solution"},n=Math.min(window.innerWidth,700),i={title:`line plot - ${a}`,width:Math.min(n,600),height:300,xaxis:{title:"x"},yaxis:{title:"Solution"},margin:{l:50,r:50,t:50,b:50}};Plotly.newPlot(o,[t],i,{responsive:!0}),console.timeEnd("plottingTime")}else if("2D"===l&&"contour"===n){let e;e=Array.isArray(r[0])?r.map((e=>e[0])):r;let t=Math.min(window.innerWidth,700),l=Math.min(...s),d=Math.max(...s),c=Math.min(...i),u=(Math.max(...i)-c)/(d-l),m=Math.min(t,600),h={title:`${n} plot - ${a}`,width:m,height:m*u,xaxis:{title:"x"},yaxis:{title:"y",scaleanchor:"x",scaleratio:1},margin:{l:50,r:50,t:50,b:50},hovermode:"closest"},f={x:s,y:i,z:e,type:"contour",line:{smoothing:.85},contours:{coloring:"heatmap",showlabels:!1},colorbar:{title:"Solution"},name:"Solution Field"};Plotly.newPlot(o,[f],h,{responsive:!0}),console.timeEnd("plottingTime")}}function re(e,t,n,o){console.time("plottingTime");const{nodesXCoordinates:s,nodesYCoordinates:i}=t.nodesCoordinates,r=e.meshConfig.meshDimension,a=V(e.meshConfig),l=new x({meshDimension:e.meshConfig.meshDimension,elementOrder:e.meshConfig.elementOrder});if("1D"===r&&"line"===n);else if("2D"===r&&"contour"===n){const r=[],d=[],c=Math.max(...s)-Math.min(...s),u=Math.max(...i)-Math.min(...i),m=50,h=Math.round(c*m),f=Math.round(u*m),p=c/(h-1),g=u/(f-1);let b=[];r[0]=Math.min(...s),d[0]=Math.min(...i);for(let e=1;e<f;e++)r[e]=r[0],d[e]=d[0]+e*g;for(let e=1;e<h;e++){const t=e*f;r[t]=r[0]+e*p,d[t]=d[0];for(let e=1;e<f;e++)r[t+e]=r[t],d[t+e]=d[t]+e*g}b=new Array(h*f).fill(null);const y=B(a),{nodeNeighbors:v,neighborCount:E}=Y(a);let C=0;for(let n=0;n<h*f;n++){if(!de(r[n],d[n],y))continue;let o=!1;for(let s=0;s<a.nop[C].length;s++){let i=a.nop[C][s]-1;for(let s=0;s<E[i];s++){let c=v[i][s];const u=ae(e,a,t,c,r[n],d[n],l);if(u.inside){C=c,b[n]=u.value,o=!0;break}}if(o)break}if(!o)for(let s=0;s<a.nop.length;s++){const i=ae(e,a,t,s,r[n],d[n],l);if(i.inside){C=s,b[n]=i.value,o=!0;break}}}let D=Math.min(window.innerWidth,700),M=u/c,F=Math.min(D,600),$=F*M,w={title:`${n} plot (interpolated) - ${e.solverConfig}`,width:F,height:$,xaxis:{title:"x"},yaxis:{title:"y",scaleanchor:"x",scaleratio:1},margin:{l:50,r:50,t:50,b:50},hovermode:"closest"},A={x:r,y:d,z:b,type:"contour",connectgaps:!1,hoverongaps:!1,line:{smoothing:.85},contours:{coloring:"heatmap",showlabels:!1},colorbar:{title:"Solution"},name:"Interpolated Solution Field"};Plotly.newPlot(o,[A],w,{responsive:!0}),console.timeEnd("plottingTime")}}function ae(e,t,n,o,s,i,r){const{nodesXCoordinates:a,nodesYCoordinates:l}=n.nodesCoordinates,d=t.nop[o].length;if(4===d){const d=R(s,i,[[a[t.nop[o][0]-1],l[t.nop[o][0]-1]],[a[t.nop[o][1]-1],l[t.nop[o][1]-1]],[a[t.nop[o][2]-1],l[t.nop[o][2]-1]],[a[t.nop[o][3]-1],l[t.nop[o][3]-1]]]);if(d.inside)return{inside:!0,value:le(e,t,n,o,d.ksi,d.eta,r)}}else if(9===d){const d=R(s,i,[[a[t.nop[o][0]-1],l[t.nop[o][0]-1]],[a[t.nop[o][2]-1],l[t.nop[o][2]-1]],[a[t.nop[o][6]-1],l[t.nop[o][6]-1]],[a[t.nop[o][8]-1],l[t.nop[o][8]-1]]]);if(d.inside)return{inside:!0,value:le(e,t,n,o,d.ksi,d.eta,r)}}return{inside:!1,value:null}}function le(e,t,n,o,s,i,r){const a=n.solutionVector,l=t.nop[o].length;let d,c=r.getBasisFunctions(s,i).basisFunction;d=Array.isArray(a[0])?a.map((e=>e[0])):a;let u=0;for(let e=0;e<l;e++)u+=d[t.nop[o][e]-1]*c[e];return u}function de(e,t,n){let o=!1;for(let s=0;s<n.length;s++){const[[i,r],[a,l]]=n[s];r>t!=l>t&&e<(a-i)*(t-r)/(l-r)+i&&(o=!o)}return o}let ce=null;async function ue(){if(ce)return ce;await import("@kitware/vtk.js/Rendering/Profiles/Geometry.js");const[{default:e},{default:t},{default:n},{default:o},{default:s},{default:i},{default:r},{default:a},{default:l},{default:d}]=await Promise.all([import("@kitware/vtk.js/Rendering/Core/Actor.js"),import("@kitware/vtk.js/Rendering/Core/ColorTransferFunction.js"),import("@kitware/vtk.js/Rendering/Core/ColorTransferFunction/ColorMaps.js"),import("@kitware/vtk.js/Common/Core/DataArray.js"),import("@kitware/vtk.js/Common/DataModel/ImageData.js"),import("@kitware/vtk.js/Filters/General/ImageMarchingSquares.js"),import("@kitware/vtk.js/Rendering/Misc/GenericRenderWindow.js"),import("@kitware/vtk.js/Rendering/Core/Mapper.js"),import("@kitware/vtk.js/Common/DataModel/PolyData.js"),import("@kitware/vtk.js/Rendering/Core/ScalarBarActor.js")]);return ce={vtkActor:e,vtkColorTransferFunction:t,vtkColorMaps:n,vtkDataArray:o,vtkImageData:s,vtkImageMarchingSquares:i,vtkGenericRenderWindow:r,vtkMapper:a,vtkPolyData:l,vtkScalarBarActor:d},ce}const me=new Map;function he(e={}){return{presetName:e.presetName??"Cool to Warm",reverse:e.reverse??!1,showScalarBar:e.showScalarBar??!0,scalarBarTitle:e.scalarBarTitle??"Solution"}}function fe(e={}){return{enabled:e.enabled??!0,numberOfContours:e.numberOfContours??12,color:e.color??[.15,.15,.15],lineWidth:e.lineWidth??1}}async function pe(e,t,n,o,s={}){console.time("plottingTime");const i=e.meshConfig.meshDimension,r=V(e.meshConfig),a=await be(e,t,r,{mode:"1D"===i&&"line"===n?"line":"surface"});await Ee(a,o,e.solverConfig,n,s),console.timeEnd("plottingTime")}async function ge(e,t,n,o,s={}){if(console.time("plottingTime"),"2D"!==e.meshConfig.meshDimension||"contour"!==n){const i=V(e.meshConfig),r=await be(e,t,i,{mode:"1D"===e.meshConfig.meshDimension&&"line"===n?"line":"surface"});return await Ee(r,o,e.solverConfig,n,s),void console.timeEnd("plottingTime")}const i=V(e.meshConfig),r=await async function(e,t,n){const o=await async function(e,t,n){const{nodesXCoordinates:o,nodesYCoordinates:s}=t.nodesCoordinates,i=new x({meshDimension:e.meshConfig.meshDimension,elementOrder:e.meshConfig.elementOrder});let r=o[0],a=o[0],l=s[0],d=s[0];for(let e=1;e<o.length;e++){const t=o[e],n=s[e];t<r&&(r=t),t>a&&(a=t),n<l&&(l=n),n>d&&(d=n)}const c=a-r,u=d-l,m=50,h=Math.max(2,Math.round(c*m)),f=Math.max(2,Math.round(u*m)),p=c/(h-1),g=u/(f-1),b=h*f,y=new Float32Array(b),v=new Float32Array(b),E=new Float32Array(b);E.fill(Number.NaN);const C=new Uint8Array(b),D=B(n),{nodeNeighbors:M,neighborCount:F}=Y(n);let $=0;for(let o=0;o<h;o++)for(let s=0;s<f;s++){const a=o*f+s,d=r+o*p,c=l+s*g;if(y[a]=d,v[a]=c,!Ne(d,c,D))continue;let u=!1;for(let o=0;o<n.nop[$].length;o++){const s=n.nop[$][o]-1;for(let o=0;o<F[s];o++){const r=M[s][o],l=Ae(e,n,t,r,d,c,i);if(l.inside){$=r,E[a]=l.value,C[a]=1,u=!0;break}}if(u)break}if(!u)for(let o=0;o<n.nop.length;o++){const s=Ae(e,n,t,o,d,c,i);if(s.inside){$=o,E[a]=s.value,C[a]=1;break}}}return{visNodesX:h,visNodesY:f,minX:r,minY:l,deltaX:p,deltaY:g,lengthX:c,lengthY:u,visNodeXCoordinates:y,visNodeYCoordinates:v,visSolution:E,insideMask:C}}(e,t,n),{visNodesX:s,visNodesY:i,minX:r,minY:a,deltaX:l,deltaY:d,visNodeXCoordinates:c,visNodeYCoordinates:u,visSolution:m,insideMask:h}=o,f=De(c,u),p=function(e,t,n){const o=[];for(let s=0;s<e-1;s++)for(let e=0;e<t-1;e++){const i=s*t+e,r=(s+1)*t+e,a=(s+1)*t+(e+1),l=s*t+(e+1);n[i]&&n[r]&&n[a]&&n[l]&&o.push(4,i,r,a,l)}return Uint32Array.from(o)}(s,i,h),g=await Ce(f,m,p,"surface");return{points:f,scalars:m,cells:p,polyData:g,mode:"surface",metadata:{meshDimension:"2D",numberOfPoints:f.length/3,numberOfCells:$e(p),interpolationGrid:{nx:s,ny:i,origin:[r,a],spacing:[l,d],imageScalars:we(m,s,i)}}}}(e,t,i);await Ee(r,o,e.solverConfig,`${n}-interpolated`,s),console.timeEnd("plottingTime")}async function be(e,t,n=null,o={}){const s=n??V(e.meshConfig),{nodesXCoordinates:i,nodesYCoordinates:r}=t.nodesCoordinates,a=Me(t.solutionVector,i.length),l=De(i,r),d=o.mode??"surface",c="line"===d?function(e){if(e<2)return new Uint32Array(0);const t=new Uint32Array(3*(e-1));let n=0;for(let o=0;o<e-1;o++)t[n++]=2,t[n++]=o,t[n++]=o+1;return t}(i.length):function(e){const t=[];for(let n=0;n<e.length;n++){const o=Fe(e[n]);t.push(o.length,...o)}return Uint32Array.from(t)}(s.nop??[]);return{points:l,scalars:a,cells:c,polyData:await Ce(l,a,c,d),mode:d,metadata:{solverConfig:e.solverConfig,meshDimension:e.meshConfig.meshDimension,numberOfPoints:l.length/3,numberOfCells:$e(c)}}}async function ye(e,t,n=null,o={}){return function(e){const{connectivity:t,offsets:n}=function(e){const t=[],n=[];let o=0,s=0;for(;o<e.length;){const i=e[o++];for(let n=0;n<i;n++)t.push(e[o++]);s+=i,n.push(s)}return{connectivity:t,offsets:n}}(e.cells),o=e.points.length/3,s="line"===e.mode,i=s?"Lines":"Polys";return['<?xml version="1.0"?>','<VTKFile type="PolyData" version="0.1" byte_order="LittleEndian">',"  <PolyData>",`    <Piece NumberOfPoints="${o}" NumberOfVerts="0" NumberOfLines="${s?n.length:0}" NumberOfStrips="0" NumberOfPolys="${s?0:n.length}">`,'      <PointData Scalars="solution">',`        <DataArray type="Float32" Name="solution" NumberOfComponents="1" format="ascii">${Array.from(e.scalars).join(" ")}</DataArray>`,"      </PointData>","      <Points>",`        <DataArray type="Float32" NumberOfComponents="3" format="ascii">${Array.from(e.points).join(" ")}</DataArray>`,"      </Points>",`      <${i}>`,`        <DataArray type="Int32" Name="connectivity" format="ascii">${t.join(" ")}</DataArray>`,`        <DataArray type="Int32" Name="offsets" format="ascii">${n.join(" ")}</DataArray>`,`      </${i}>`,"    </Piece>","  </PolyData>","</VTKFile>"].join("\n")}(await be(e,t,n,o))}function ve(e){const{nodesXCoordinates:t,nodesYCoordinates:n}=e.nodesCoordinates,o=Me(e.solutionVector,t.length),s=new Float32Array(3*t.length);for(let e=0;e<t.length;e++){const i=3*e;s[i]=Number(t[e])||0,s[i+1]=Number(n[e])||0,s[i+2]=o[e]}return{features:s,featuresShape:[t.length,3],labels:o,points:De(t,n)}}async function Ee(e,t,n,o,s={}){if("undefined"==typeof document)return void i("vtk.js visualization requires a browser environment");const{vtkActor:r,vtkColorTransferFunction:a,vtkColorMaps:l,vtkGenericRenderWindow:d,vtkMapper:c,vtkScalarBarActor:u}=await ue(),m=document.getElementById(t);if(!m)return void i(`Could not find plot container with id: ${t}`);me.has(t)&&(me.get(t).delete(),me.delete(t)),m.innerHTML="",m.style.position="relative",m.style.width=m.style.width||"100%",m.style.height=m.style.height||"420px";const h=d.newInstance({background:[1,1,1]});h.setContainer(m),h.resize();const f=h.getRenderer(),p=h.getRenderWindow(),g=c.newInstance();g.setInputData(e.polyData),g.setScalarModeToUsePointData(),g.setColorModeToMapScalars(),g.setScalarVisibility(!0);const b=function(e){if(!e?.length)return[0,1];let t=Number.POSITIVE_INFINITY,n=Number.NEGATIVE_INFINITY;for(let o=0;o<e.length;o++){const s=e[o];Number.isFinite(s)&&(s<t&&(t=s),s>n&&(n=s))}return Number.isFinite(t)&&Number.isFinite(n)?t===n?[t-1,n+1]:[t,n]:[0,1]}(e.scalars);g.setScalarRange(b[0],b[1]);const y=s.colorScale??he({}),v=a.newInstance(),E=function(e,t){if(!t||!e?.RGBPoints?.length)return e;const n=e.RGBPoints,o=n[0],s=n[n.length-4],i=[];for(let e=n.length-4;e>=0;e-=4)i.push(s-(n[e]-o),n[e+1],n[e+2],n[e+3]);return{...e,RGBPoints:i}}(l.getPresetByName(y.presetName)??l.getPresetByName("Cool to Warm"),y.reverse);v.applyColorMap(E),v.setMappingRange(b[0],b[1]),v.updateRange(),g.setLookupTable(v);const C=r.newInstance();if(C.setMapper(g),"line"===e.mode&&C.getProperty().setLineWidth(3),f.addActor(C),y.showScalarBar){const e=u.newInstance({drawNanAnnotation:!1,generateTicks:e=>{const t=e.getLastTickBounds();if(!t||t.length<2)return;const[n,o]=t,s=(o-n)/4,i=Array.from({length:5},((e,t)=>n+t*s));e.setTicks(i),e.setTickStrings(i.map((e=>{const t=Math.abs(e);return 0===t?"0":t>=.01&&t<1e4?parseFloat(e.toPrecision(4)).toString():e.toExponential(2)})))}});e.setTickTextStyle({fontColor:"black"}),e.setAxisTextStyle({fontColor:"black"}),e.setAxisLabel(y.scalarBarTitle),e.setScalarsToColors(v),f.addActor2D(e)}const D=fe(s.contourLines??{enabled:!1});D.enabled&&"line"!==e.mode&&await async function(e,t,n,o){const s=t.metadata?.interpolationGrid;if(!s)return;const{vtkActor:i,vtkDataArray:r,vtkImageData:a,vtkImageMarchingSquares:l,vtkMapper:d}=await ue(),c=a.newInstance();c.setDimensions(s.nx,s.ny,1),c.setOrigin(s.origin[0],s.origin[1],0),c.setSpacing(s.spacing[0],s.spacing[1],1);const u=r.newInstance({name:"solution",numberOfComponents:1,values:s.imageScalars});c.getPointData().setScalars(u);const m=l.newInstance({slicingMode:2,slice:0,mergePoints:!0});m.setInputData(c);const h=Math.max(2,o.numberOfContours),f=n[0],p=n[1],g=(p-f)/(h-1),b=[];for(let e=0;e<h;e++)b.push(f+e*g);m.setContourValues(b),m.update();const y=d.newInstance();y.setInputData(m.getOutputData()),y.setScalarVisibility(!1);const v=i.newInstance();v.setMapper(y),v.getProperty().setColor(...o.color),v.getProperty().setLineWidth(o.lineWidth),e.addActor(v)}(f,e,b,D),f.resetCamera(),p.render(),me.set(t,h),m.title=`${o} plot - ${n}`}async function Ce(e,t,n,o="surface"){const{vtkPolyData:s,vtkDataArray:i}=await ue(),r=s.newInstance();r.getPoints().setData(e,3),n.length>0&&("line"===o?r.getLines().setData(n):r.getPolys().setData(n));const a=i.newInstance({name:"solution",numberOfComponents:1,values:t});return r.getPointData().setScalars(a),r}function De(e,t){const n=new Float32Array(3*e.length);for(let o=0;o<e.length;o++){const s=3*o;n[s]=Number(e[o])||0,n[s+1]=Number(t?.[o])||0,n[s+2]=0}return n}function Me(e,t){const n=new Float32Array(t);for(let o=0;o<t;o++){const t=e?.[o];n[o]=Number(Array.isArray(t)?t[0]:t)||0}return n}function Fe(e){const t=e.map((e=>e-1)),n=t.length;return 2===n||3===n?t:4===n?[t[0],t[2],t[3],t[1]]:6===n?[t[0],t[2],t[5]]:8===n?[t[0],t[4],t[6],t[2]]:9===n?[t[0],t[6],t[8],t[2]]:t.slice(0,Math.min(4,t.length))}function $e(e){let t=0,n=0;for(;n<e.length;)n+=e[n]+1,t++;return t}function we(e,t,n){const o=new Float32Array(t*n);for(let s=0;s<n;s++)for(let i=0;i<t;i++)o[i+s*t]=e[i*n+s];return o}function Ae(e,t,n,o,s,i,r){const{nodesXCoordinates:a,nodesYCoordinates:l}=n.nodesCoordinates,d=t.nop[o].length;if(4===d){const d=R(s,i,[[a[t.nop[o][0]-1],l[t.nop[o][0]-1]],[a[t.nop[o][1]-1],l[t.nop[o][1]-1]],[a[t.nop[o][2]-1],l[t.nop[o][2]-1]],[a[t.nop[o][3]-1],l[t.nop[o][3]-1]]]);if(d.inside)return{inside:!0,value:xe(e,t,n,o,d.ksi,d.eta,r)}}else if(9===d){const d=R(s,i,[[a[t.nop[o][0]-1],l[t.nop[o][0]-1]],[a[t.nop[o][2]-1],l[t.nop[o][2]-1]],[a[t.nop[o][6]-1],l[t.nop[o][6]-1]],[a[t.nop[o][8]-1],l[t.nop[o][8]-1]]]);if(d.inside)return{inside:!0,value:xe(e,t,n,o,d.ksi,d.eta,r)}}return{inside:!1,value:null}}function xe(e,t,n,o,s,i,r){const a=n.solutionVector,l=t.nop[o].length,d=r.getBasisFunctions(s,i).basisFunction,c=Array.isArray(a[0])?a.map((e=>e[0])):a;let u=0;for(let e=0;e<l;e++)u+=c[t.nop[o][e]-1]*d[e];return u}function Ne(e,t,n){let o=!1;for(let s=0;s<n.length;s++){const[[i,r],[a,l]]=n[s];r>t!=l>t&&e<(a-i)*(t-r)/(l-r)+i&&(o=!o)}return o}class ke{constructor(){this.worker=null,this.feaWorker=null,this.isReady=!1,this._initWorker()}async _initWorker(){try{this.worker=new Worker(new URL("./wrapperScript.js",import.meta.url),{type:"module"}),this.worker.onerror=e=>{console.error("FEAScriptWorker: Worker error:",e)};const e=p(this.worker);this.feaWorker=await new e,this.isReady=!0}catch(e){throw console.error("Failed to initialize worker",e),e}}async _ensureReady(){return this.isReady?Promise.resolve():new Promise(((e,t)=>{let n=0;const o=()=>{n++,this.isReady?e():n>=50?t(new Error("Timeout waiting for worker to be ready")):setTimeout(o,1e3)};o()}))}async setModelConfig(e){return await this._ensureReady(),this.feaWorker.setModelConfig(e)}async setMeshConfig(e){return await this._ensureReady(),this.feaWorker.setMeshConfig(e)}async addBoundaryCondition(e,t){return await this._ensureReady(),this.feaWorker.addBoundaryCondition(e,t)}async setSolverMethod(e){return await this._ensureReady(),this.feaWorker.setSolverMethod(e)}async solve(){await this._ensureReady(),performance.now();const e=await this.feaWorker.solve();return performance.now(),e}async getModelInfo(){return await this._ensureReady(),this.feaWorker.getModelInfo()}async ping(){return await this._ensureReady(),this.feaWorker.ping()}terminate(){this.worker&&(this.worker.terminate(),this.worker=null,this.feaWorker=null,this.isReady=!1)}}const Pe="0.3.0 (RC)";export{oe as FEAScriptModel,ke as FEAScriptWorker,he as createColorScale,fe as createContourLineOptions,se as importGmshQuadTri,n as logSystem,re as plotInterpolatedSolution,ge as plotInterpolatedSolutionVtk,ie as plotSolution,pe as plotSolutionVtk,Pe as printVersion,ve as transformSolverOutputToMLBuffers,ye as transformSolverOutputToVTP,be as transformSolverOutputToVtkData};
//# sourceMappingURL=feascript.esm.js.map