1. A method of modeling thermal problems using a non-dimensional finite element method, comprising a computer executing the steps of: (a) establishing a steady state heat conduction equation in three dimensions for an orthotropic material as k x ⁢ ∂ 2 ⁢ T ∂ x 2 + k y ⁢ ∂ 2 ⁢ T ∂ y 2 + k z ⁢ ∂ 2 ⁢ T ∂ z 2 + G = 0 , wherein the orthotropic material has three surfaces having surface areas S 1 , S 2 and S 3 , respectively, wherein T=T b holds for S 1 , k x ⁢ ∂ T ∂ x ⁢ l ^ + k y ⁢ ∂ T ∂ y ⁢ m ^ + k z ⁢ ∂ T ∂ z ⁢ n ^ + q * = 0 holds for S 2 , and k x ⁢ ∂ T ∂ x ⁢ l ^ + k y ⁢ ∂ T ∂ y ⁢ m ^ + k z ⁢ ∂ T ∂ z ⁢ n ^ + h ⁡ ( T - T ∞ ) = 0 holds for S 3 , where k x , k y , and k z are thermal conductivities in Cartesian x, y and z-directions, respectively, G is a rate of internal heat generation per unit volume, T is temperature, T b and T ∞ are specified and ambient temperatures, respectively, q* is specified heat flux, h is a coefficient of convective heat transfer, and {circumflex over (l)}, {circumflex over (m)} and {circumflex over (n)} are surface normals, respectively corresponding to the Cartesian x, y and z-directions; (b) establishing a scalar functional quantity Π=U+Ω G +Ω q +Ω h , wherein U represents internal energy, Ω G represents internal heat generation, Ω h represents heat convection, and Ω q represents heat conduction, and: U = 1 2 ⁢ ∫ ∫ ∫ V ⁡ [ k x ⁡ ( ∂ T ∂ x ) 2 + k y ⁡ ( ∂ T ∂ y ) 2 + k z ⁡ ( ∂ T ∂ z ) 2 ] ⁢ ⅆ V , ⁢ Ω G = - ∫ ∫ ∫ V ⁢ GT ⁢ ⅆ V , ⁢ Ω q = - ∫ ∫ ∫ S ⁢ ⁢ 2 ⁢ q * ⁢ T ⁢ ⅆ S , and Ω h = 1 2 ⁢ ∫ ∫ S ⁢ ⁢ 3 ⁢ h ⁡ ( T - T ∞ ) 2 ⁢ ⅆ S , wherein V represents volume; (c) establishing a set of non-dimensional parameters θ = T T ∞ , x _ = x L x , y _ = y L y , and ⁢ ⁢ z _ = z L z , where θ is non-dimensional temperature, x , y and z are non-dimensional spatial variables, corresponding to a domain to be discretized and L x , L y and L Z are maximum dimensions of a domain to be discretized along the x, y and z-directions, respectively; (d) expressing the scalar functional quantity for no internal heat generation in matrix form as: Π hT ∞ 2 = 1 2 ⁢ ∫ ∫ ∫ V ⁢ { g _ } T ⁡ [ L _ ] ⁢ { g _ } ⁢ ⅆ V - ∫ ∫ S ⁢ ⁢ 2 ⁢ { θ } T ⁡ [ N _ ] ⁢ T ⁢ q _ * ⁢ ⅆ S + 1 2 ⁢ ∫ ∫ S ⁢ ⁢ 3 ⁢ ( { θ } T ⁡ [ N _ ] T - 1 ) 2 ⁢ ⅆ S , where [ N ] is a shape function matrix, { g } is a non-dimensional temperature gradient, {θ} is a non-dimensional temperature vector, [ D ] is a non-dimensional material property matrix, and [ L ] is a geometric dimensions ratio matrix, wherein { g _ } = [ ∂ θ ∂ x _ ∂ θ ∂ y _ ∂ θ ∂ z _ ] T ⁢ ⁢ and { g _ } = [ B _ ] ⁢ { θ } , where ⁢ [ B _ ] = [ ∂ [ N _ ] ∂ x _ ∂ [ N _ ] ∂ y _ ∂ [ N _ ] ∂ z _ ] T , and [ D _ ] = [ 1 Bi x 0 0 0 1 Bi y 0 0 0 1 Bi z ] , [ L _ ] = [ 1 0 0 0 L x L y 0 0 0 L x L z ] , where Bi x , Bi y and Bi z are Biot numbers given by Bi x = hL x k x , Bi y = hL y k y ⁢ ⁢ and ⁢ ⁢ Bi z = hL z k z and non-dimensional specified heat flux q * is given by q _ * = q * hT ∞ ; (e) minimizing a dimensionless scalar functional quantity given as Π=Π/hT ∞ 2 L y L z with respect to {θ} to yield: [ ∫ ∫ ∫ V _ ⁡ [ B _ ] T ⁡ [ L _ ] ⁡ [ D _ ] ⁡ [ B _ ] ⁢ ⅆ V _ + ∫ ∫ S _ ⁢ ⁢ 3 ⁡ [ [ N _ ] T ⁡ [ N _ ] ] ⁢ ⅆ S _ ] ⁢ { θ } = ∫ ∫ S _ ⁢ ⁢ 2 ⁡ [ N _ ] T ⁢ q _ * ⁢ ⅆ S _ + ∫ ∫ S _ ⁢ ⁢ 3 ⁡ [ N _ ] T ⁢ ⅆ S _ ; (f) calculating an element stiffness matrix as [ k ]=[ k q ]+[ k h ], wherein [ k _ q ] = ∫ ∫ ∫ V _ ⁡ [ B _ ] T ⁡ [ L _ ] ⁡ [ D _ ] ⁡ [ B _ ] ⁢ ⅆ V _ ⁢ ⁢ and ⁢ [ k _ h ] = ∫ ∫ S _ ⁢ ⁢ 3 ⁢ ⌊ [ N _ ] T ⁡ [ N _ ] ⌋ ⁢ ⅆ S _ ; and (g) calculating an element load vector as { f }={ f q }+{ f h }, wherein: [ f _ q ] = ∫ ∫ S _ ⁢ ⁢ 2 ⁡ [ N _ ] T ⁢ q _ * ⁢ ⅆ S _ ⁢ ⁢ and ⁢ [ f _ h ] = ∫ ∫ S _ ⁢ ⁢ 3 ⁡ [ N _ ] ⁢ ⅆ S _ .