Anisotropic Mesh Adaptation Based on Error Estimation for 3-D Finite-Element Simulation of Electromagnetic Material Processing

This work aims to develop an adaptive remeshing procedure for the finite-element method (FEM) on electromagnetic computations. A thorough comparison of metric computation strategies is carried out as it constitutes a cornerstone of this development. This procedure will focus on the mesh size adaptation to distribute the error uniformly over a computational domain, to achieve a user-prescribed solution accuracy. Also, it shall enable dealing with complex geometries for electromagnetic-coupled material processing applications. For this purpose, a quasi-steady-state approximation of Maxwell’s equations in a time-domain formalism is considered. The automatic remeshing procedure is based on the following key steps: An a posteriori error estimator to pinpoint the critical areas needing refinement or allowing coarsening. An anisotropic metric approximation. Both steps use a global field recovery algorithm to enable robust gradient computation. Finally, several 3-D test cases are presented.


I. INTRODUCTION
E LECTROMAGNETIC material processing often involves parts with complex geometries and materials with nonlinear magnetic behavior. Efficient modeling of these processes can often be highly demanding in terms of CPU time and resources. Anisotropic mesh adaptation allows capturing the behavior of complex physical phenomena as well as improving the accuracy of the numerical solution. Appropriate control of the algorithm ensures that an adequate number of degrees of freedom are maintained to achieve the desired accuracy. This work aims to construct a metric tensor to carry on an anisotropic mesh adaptation procedure dealing with the electromagnetic phenomenon. It will be applied to a fully immersed computational model that we have developed, thus enabling all conductive and nonconductive materials as well as the surrounding domain (air) to be remeshed.

A. Electromagnetic-Coupled Processes and Modeling
Electromagnetic processing of materials (EPM) has gained popularity in recent decades, becoming a large family of industrial technologies due to its many advantages. These processes are often more energy-efficient than a purely thermomechanical processing method. They offer many advantages, including rapid heat-up and production startup speed, minimal material waste, and high-energy transfer efficiency. Among the most popular, we can mention induction heating (IH) [1], [2], which is nowadays one of the most widely used, allows heat The continuous improvement of computational processing units (CPUs/GPUs) and the availability of high-performance computing (HPC) infrastructures allow computational models to become increasingly realistic and accurate, which in turn means increased numerical complexity. To reduce the time from concept to production by improving the engineering design phase, one of the main challenges and opportunities for improving numerical methods is to automatize the generation and adaptation of meshes. In the FEM, the computational domain is divided into a finite number of elements on which the system of partial differential equations describing the physical phenomena is discretized. To achieve an accurate numerical solution, elements must have a minimum geometrical quality to limit the influence of round-off errors intrinsic to finite precision arithmetic. Moreover, a very fine discretization results in a large number of elements leading to large memory and CPU time requirements. To overcome this situation, mesh adaptation methods have been developed. This involves locating the areas where the numerical solution is not sufficiently accurate and proposing local refinements as well as de-refinements if the resources can be located elsewhere in the domain.

B. Metric-Based Anisotropic Mesh Adaptation-Basic Principles
The link between mesh generation and a suitable numerical solution has been a strong research topic for several years in the engineering community. Enabling a sufficient accuracy of the numerical solution with a reduced number of degrees of freedom in the mesh is the main objective to reduce the computational burden.
Anisotropic mesh adaptation is a useful tool widely used in mesh generation to achieve the best possible solution. It allows adaptation of the size, shape, and orientation of an element according to a specific quantity that can improve the accuracy of the solution, while improving the computational performance. Anisotropic mesh adaptation aims to create an optimal mesh of the computational domain in the sense that the error of the solution of the considered problem is uniformly distributed over the entire domain.
The main goal of this approach is the correct calculation of the metric, which is usually based on some error estimate. It is worth mentioning that error estimators based on a local error problem can be inaccurate in producing anisotropic meshes [13]. In fact, they usually do not contain enough information about the direction of the solution, which is global in nature, and their accuracy and efficiency are sensitive to the aspect ratio of the elements, which can be large for anisotropic meshes. Therefore, it is desirable to define an error estimate using a globally defined error problem. Some methods have been published on the use of a posteriori error estimators to control anisotropic mesh adaptation. Huang et al. [14] developed and compared a hierarchical error estimator with a recovery-based method. Cao et al. [15] calculated some scalar monitoring functions to use adaptive mesh motion from interpolation errors and a posteriori error estimates. Apel et al. [16] used the gradients of some a posteriori error estimates for calculating anisotropic meshes.
Agouzal et al. [17] calculated a metric tensor from edge-based error estimates.
One approach to constructing the metric tensor field, for the generation of adaptive mesh, consists of computing the Hessian of the exact solution of a scalar field. Since finding the exact solution is, in most cases, impossible, numerical methods must be used to construct it from the finite-element solution. One of the most commonly used methods is based on recovery strategies. In fact, these strategies are advantageous since they are not tied to a particular numerical scheme or a specific application and can be applied to a wide range of PDEs, providing high-dimensional information about the solution [18], [19]. Most of these recovery techniques, which typically use a discrete least squares or integral formulation, compute nodal values that are convex combinations of the gradient over neighboring elements. Other recovery techniques include simple weighted average methods [20] and minimization problems [21]. Lipnikov and Vassilevski [22] compared these methods by comparing three of them: a variational method [23] and two projection methods [24], [25]. They concluded that the variational method is less expensive and offers the most accurate interpolation solution as well as the most robust behavior over a larger scale of norms than the other methods. Another variational method for Hessian recovery was introduced by Agouzal and Vassilevski [26].
In electromagnetic modeling, adaptive mesh has been little studied so far, we can mention the work of Grosges et al. [28] and [29] in 2-D and 3-D for simple geometries in the calculation of high electromagnetic field gradients [27], [28], [29]. In this work, we develop an adaptive remeshing procedure for FEMs, based on Hessian or Jacobian recovery strategy, which allows us to deal with complex geometries for electromagnetically coupled material processing applications.
The structure of the article is as follows. Section II presents the physical phenomenon and its intrinsic anisotropy. Section III presents two strategies for calculating the metric field. Section IV compares two different ways of symmetrizing the tensor field. Section V presents numerical results for both metrics in different electromagnetic processes. Conclusions are given in Section VI.
II. ELECTROMAGNETIC PHENOMENA In electromagnetic materials processing, various phenomena can produce a high degree of anisotropy. Depending on the processing conditions, abrupt changes in material properties can occur. The skin effect and the Curie temperature are the main sources of anisotropy in electromagnetic phenomena and therefore have a major influence on these processes. Therefore, the main objective of this work is to develop a relatively general approach that is independent of these phenomena and can capture them. We start by introducing the equations describing the electromagnetic phenomenon.

A. Generalization for Electromagnetic Modeling
The physical phenomena of magnetic and electric wave propagation are described by Maxwell's equations where: E: electric field intensity; D: electric flux density; H : magnetic field intensity; B: magnetic flux density; J : electric current density; and ρ: electric charge density.
In this system of equations, relations (1) and (2) express the coupling between the electric ⃗ E, ⃗ D and magnetic ⃗ B, ⃗ H fields. Most electromagnetic material processes usually take place at low/medium frequencies (<1 MHz), thus enabling to neglect the current displacements density term ∂ t ⃗ D (quasisteady-state approximation). Then, for our purpose, Maxwell's equations can be rewritten as Current density conservation is expressed as The previous equations are not complete, and the Maxwell equations involve only the fields and their sources. As in solid mechanics, constitutive relations, modeling the electromagnetic properties of the different media, must be added. The constitutive relations describing the macroscopic properties of the medium being considered are: Magnetic permeability (µ) Electrical conductivity (Ohm's law) The proportionality parameter σ is called the electrical conductivity.
The electromagnetic problem of (5)-(7) is solved in a A −φ formulation on time-domain discretization by the FEM. For further details, the reader is referred to [30] and [31].

B. Skin Effect
Most electromagnetic-coupled processes deal with timevarying electrical loadings (harmonic or pulsed). As a time-dependent current runs through the coils, a time-varying electromagnetic field is induced in the surrounding domain. This electromagnetic field generates eddy currents in the conductive workpiece. These currents dissipate heat through the Joule effect and produce Lorentz forces, thus allowing the workpiece to be heated and deformed.
The distribution of eddy currents circulation within the workpiece is, however, not uniform. The highest density is reached on the surface and decreases when going toward the part's core (Fig. 1). This phenomenon is called "the skin depth In dc current flows homogeneously on the cross section. In ac the current alternates not only in time, but also its peakdensity location on the cross section. effect" and for perfectly harmonic loading, it is described by the following expression: where f represents the signal's frequency, ρ the material's electrical resistivity, and µ the magnetic permeability (index 0 stands for void absolute permeability and r the material's relative value). It is very important in an IH system to control these parameters, both for surface treatments as well as for mass heating. Low frequencies enable greater heating depths, which are useful to heat massive parts; higher frequencies induce currents that are mainly concentrated on the surface. Thanks to this specificity, this process is commonly used in surface heat treatment processes.

C. Curie Temperature
Another source of anisotropy is the Curie point, which can drastically change material properties at a given temperature. The magnetic properties of materials come from different intrinsic magnetic moment structures, which are directly dependent on temperature. The magnetic moment is basically the magnetic strength and orientation of a magnet or any other object that produces a magnetic field. Depending on their arrangement, materials can be classified as ferromagnetic and paramagnetic [ Fig. 2(b)]. Ferromagnetic materials possess equal parallel magnetic moments, produced by the alignment of the magnetic dipoles within the material and are ordered on the same magnitude in the absence of an applied magnetic field. On the other hand, the magnetic moments of a paramagnetic material are disordered in the absence of an applied magnetic field and ordered when one is applied.
Ferromagnetic materials have an inherent magnetization M, from which relation (9) is rewritten as where µ 0 is the void's magnetic permeability and χ is the magnetic susceptibility. From this latter value, the relative magnetic permeability µ r can be defined as For these materials, the relationship between the magnetic flux density B and applied magnetic field H is highly nonlinear as can be seen in Fig. 2(a) (M s is magnetization saturation), which is not the case for paramagnetic materials.
When materials are heated, high-energy electrons are produced that can disrupt the order and destroy the alignment between the dipoles, making them chaotic. This occurs when the material reaches the Curie temperature (T C ), which is the critical point at which the intrinsic magnetic moments of a material change direction. For example, the ordered magnetic moments of ferromagnets change and become disordered as paramagnetic at this temperature. Thus, losing their magnetization as shown in Fig. 3.
It should be noted that our study focuses only on the macroscopic scale of the physical phenomenon. We use this analogy on the microscopic scale to explain the source of the change in the magnetic properties of materials. In the simulation of electromagnetic phenomena, this change can be easily identified as the variation of the concentration of magnetic field lines underneath the surface. While paramagnetic materials have no effect on the flow of the magnetic field, ferromagnetic materials flatten the electromagnetic field lines; thus, depending on the value of the frequency, these field lines become increasingly parallel and close to the workpiece surface. Fig. 4 shows the magnetic field lines generated by a centrally placed inductor placed in the red square and its distribution all over the domain; the vertical red line represents the outer surface of a workpiece close to the inductor. Fig. 4(a) and (b) shows how the magnetic field lines get distorted when changing from a nonmagnetic to a ferromagnetic material.

III. ANISOTROPIC MESH ADAPTATION
Many mesh adaptation strategies lead to the construction of anisotropic meshes: metric-based mesh adaptation is the most popular method used to drive an anisotropic mesh generation procedure. Its main feature is the construction of a metric tensor that contains the necessary information to modify the mesh in an anisotropic way. To generate anisotropic meshes, it must be possible to determine at each point of the domain the desired sizes and directions of the final mesh elements. To do this, the tensor field at each point of the  domain is required to evaluate the lengths and angles in the mesh. An anisotropic mesh adaptation based on metric fields requires the construction of several quantities, such as an error estimator to control the mesh size and normalize the metric and on top of it, a quantity describing the deformation or distortion of the non-Euclidean manifold due to the action of the existing physical fields.

A. Error Estimator
The error in a numerical model is obtained by calculating the difference between the exact solution and a numerical one, which can be written as follows: where u is the exact solution and u h is the one resulting from the finite-element analysis. Error estimators are used to control the anisotropic mesh adaptation procedures. They allow a local modification of the mesh and govern the mesh to achieve a prescribed level of error in the L 2 norm.
Two main approaches to error estimation can be found in the literature: a priori and a posteriori. A priori error estimation requires prior knowledge of the solution before performing the actual calculation. Among the techniques, one can mention: the extrapolation of solutions obtained on sequences of progressively finer meshes (element size variation) or on sequences of meshes with shape functions of increasing order (degrees of interpolation functions), followed by a comparison of the solutions to obtain an indication of the error. A posteriori error estimators are developed during postprocessing of the physical fields of interest. Interest in a posteriori error estimation for FEMs in boundary value problems began with the pioneering work of Babuška and Rheinboldt [32]. A posteriori error estimation techniques were developed to approximate the error in energy or an energy norm in each finite element. This idea formed the basis of adaptive meshing procedures designed to control and minimize error. In this work, a gradient recovery method is used to compute the error.
1) Gradient Recovery Methods: Recovery-based techniques are based on the determination of a smoothed solution. They are based on the concept that a better approximation of the exact solution is achieved with an appropriate mesh. In most cases, the exact solution is unknown, but an attempt is made to construct a higher order than the one calculated by the finite-element solution that counts as a good approximation of the exact solution. Several strategies for constructing the enriched or smoothed solution can be found in the literature; again, it is important when constructing the smoothed solution to take into account robustness and the ability to cope with geometrically complex configurations, such as those found in industrial processes, as well as the use of fast and efficient methods. The most popular are undoubtedly the method of Zhu and Zienkiewicz [25] and Galerkin recovery or residual minimization recovery [31], where the gradient of the finite-element solution is recovered. Then the accuracy of the smoothed gradient is supposed to be higher than the one obtained by the finite-element approximation. Therefore, the evaluation of the difference between the smoothed and the direct gradient is used to compute the error estimator.
where ∇u h is the gradient calculated from the finite-element solution and G(u h ) is the enriched or smoothed gradient.
2) Error Estimator Calculation: The strategy used here is the Galerkin recovery or residual minimization recovery; it consists in building arbitrary p-order polynomials. The interpolation order p of this new polynomial should be at least one order higher than the original interpolation degree of the field at the Gauss points. In this procedure, the order of interpolation is constrained by the interpolation degree of the FE mesh used to represent the original field since it is built using the same mesh as the one required for solving the original PDE. Let ⃗ X P 0 be the field obtained from the finite-element analysis and ⃗ X P 1 is the recovered field. The method consists in solving a global minimization (or recovery) problem given by For all ⃗ ψ, edges basis functions, belonging to the Sobolev space chosen for the original PDE (H (curl) for the electromagnetic fields). This formulation leads to a linear system. After resolution, the error is obtained as (for more details, refer to [31]) where x G represents the Gauss point coordinates of each element (e). Then, the error normalization is given by where ϵ L 2 ( e ) and ϵ L 2 ( ) are the local and global errors, respectively. It is worth pointing out that local defined error estimates generally do not take into consideration the directional effect of the anisotropy of the solution. However, a globally defined error estimator, such as the one defined here, contains additional information about the solution, which allows the generation of an anisotropic mesh adaptation approach [14]. Once the error is computed as a field distribution on the current mesh, the next step consists in computing the metric tensor.
B. Metric-Based Methods 1) Metric Tensor: Anisotropic mesh generation is based on the definition at each node of a discrete metric tensor, which allows controlling the size, shape, and orientation of the mesh elements. Thus, the metric tensor is defined as a symmetric positive-definite square matrix. In other words, a metric is used to measure space. Therefore, the scalar product of two vectors x and y in accordance with a metric M can be expressed as follows: With the above expression, it is easy to define the associated non-Euclidean norm As a result of symmetry and positive definiteness, the metric tensor M can be decomposed as follows: with ⃗ e 1 , ⃗ e 2 , and ⃗ e 3 being the eigenvectors and λ 1 , λ 2 , and λ 3 the eigenvalues of M. It must be noted that M represents a matrix transformation of rotation and distortion.
2) Evaluation of First-and Second-Order Derivatives: A Galerkin recovery or residual minimization recovery procedure is adopted to calculate the metric field. The Hessian is a square matrix of second-order partial derivatives of a scalar-valued function that is commonly used to construct the metric tensor. This matrix measures the rate of change of a scalar field in its spatial vicinity. It is also a description of the local curvature of a multivariable function. It will, therefore, provide the appropriate information to adapt the mesh to a predetermined field. In this work, in addition to a metric tensor computed from a Hessian field, which generally uses a scalar field, computed from any field norm, a second method is proposed, in which the metric is computed as the Jacobian of a vector field. Both procedures are described as follows. a) Hessian from a scalar field Hess(∥x∥): Let X P 0 be a scalar continuous field, which is computed by recovery from the error estimator or any field norm ∥X ∥. A recovered field X P 1 at the nodes is calculated by a global minimization problem as ⟨ϕ, X P 1 , X P 0 ⟩ = 0 (27) where ϕ is the nodal basis function. Then, a discontinuous gradient vector field ⃗ ∇(X ) P 0 is calculated at the Gauss points by following the finite-element definition To construct the Hessian field, the above procedure is repeated. First by computing a continuous gradient field by recovery ⃗ G P 1 at the nodes from ⃗ ∇(X ) P 0 Then, using (28) to compute the discontinuous Hessian tensor field Hess P 0 at the Gauss points Finally, the recovery method is used to compute the continuous Hessian tensor field Hess P 1 at each node of the mesh ⟨ϕ, Hess P 1 − Hess P 0 ⟩ = 0.
b) Jacobian of a vector field J(x): To calculate the Jacobian "J " of a vector "x," the procedure applied is exactly as explained to compute the Hessian field. However, the gradient calculation is applied once. Since any gradient of a vector field is directly a matrix. The transpose of this matrix is taken as the Jacobian of any vector field, the Jacobian matrix can be thought of as describing the amount of stretching and rotation that a function imposes locally, and therefore this matrix is also used to compute the metric tensor.
Let ⃗ X P 0 be an elementwise vector field. Then a recovered continuous vector field ⃗ X P 1 , at the element nodes, can be computed by recovery as follows: Then, a discontinuous gradient field ⃗ ∇ ⃗ X P 0 at the Gauss points from ⃗ X P 1 is calculated by following the finite-element definition with ⃗ ∇ ⃗ X P 0 being directly a tensor field. Finally, the recovery procedure is again used to calculate the continuous Jacobian tensor field J P 1 from ⃗ ∇ ⃗ X P 0 at each node of the mesh.
3) Metric Computation: The resulting tensors obtained in Subsection III-B2 are generally not positive definite. For this reason, symmetry and positive definiteness must be enforced. Two methods of symmetrization have been tested. The first one consists of taking the average between the tensor and its transpose (since we want to symmetrize tensors, the letter T will be used to represent them) While this naïve approach is rather straightforward to implement and computationally cheap, it does not guarantee positive definiteness. Thus, a second approach to enforce both requirements can be built using the conjugate product a) Metric from the Hessian of a scalar field Hess(∥x∥): The tensor calculated in (36) leads to a square tensor. Meaning that a linearization is required, the eigenvectors and eigenvalues are computed from this resulting tensor R, = eig(Hess(∥x∥)). (37) Then, the metric M is calculated from the eigenvalues and eigenvectors as follows: On the other hand, the metric symmetrization calculated with (35) needs no further modification. In both cases, it leads to a metric proportional to the inverse of the square of the mesh size b) Metric from the Jacobian of a vector field J (x): In this case, the symmetrizing procedure defined in (36) allows us to properly approximate the Jacobian to a square of the first derivatives. The eigenvectors and eigenvalues are computed from the Jacobian by Then, to compute the metric M, no linearization is required Thus, as in the case of Hess(∥x∥), it also leads to the same proportionality between the metric and the mesh size For the Jacobian procedure, the symmetrization defined in (35) does not seem to fit with this metric calculation adequately, since it would not lead to the proportionality shown in (42) and is therefore not used in this case.

4) Metric Normalization:
In the metric-based anisotropic mesh adaptation, the eigenvalues of the metric tensor enable relating the required mesh size (h) along the eigenvector's directions with the error estimator, by using the proportionality ratio To avoid nonusable metric specifications outputs such as infinite mesh sizes in an element, caused, for instance, by highly directional fields, the eigenvalues are rescaled by applying an error tolerance (ϵ 0 ) computed from (20) and bounding the minimal (h min ) and maximal (h max ) mesh size values.λ In the case of anisotropy loss due to normalization, the ratio between the ratios between eigenvalues is stored as Then, after the normalization defined by (44). The eigenvalues are checked to see if their values are the same λ 1 =λ 2 =λ 3 ; in this case, the eigenvalues are readjusted according to (45) and (46). It allows for maintaining the anisotropy of the solution as well as having a reasonable mesh size.
Finally, the metric is deduced as For the sake of simplicity, the double-bar notation for second-order tensors will not be written for the final metric M.

C. Algorithms Proposed
The procedures used to compute the metric tensor are summarized in the next algorithms (see Tables I and II).

IV. ANISOTROPY CAPTURING COMPARISON
The metric calculated from a Hessian field (M(∥x∥) and the metric calculated from a Jacobian matrix M(x) of a vector field x are studied in this section.
First, the symmetrization procedure is tested and then its ability to capture the anisotropy is proved. To do this, a simple induction case is proposed, which is shown in Fig. 5, where the workpiece and the inductor are represented as thin plates.

A. Symmetrization Procedure
To begin, only the metric M(∥x∥) is used to compare both symmetrization procedures. Since a metric tensor can be represented as an ellipsoid in 3-D space, where the principal axes are given by their eigenvectors and their length by the eigenvalues [33]. Knowing this, and applying a current density as shown in Fig. 5, we should expect well-ordered ellipsoids in the same direction. Fig. 6 shows the ellipsoids of the metric tensor M(∥x∥) built with the magnetic field ( ⃗ H ) for both symmetrization procedure defined in (35) and (36).
For each one of the symmetrization procedures T 1 and T 2 displayed in Fig. 6(a) and (b), respectively, the behavior found is quite different. It can be seen that, although the T 1 procedure predicts a globally consistent anisotropy behavior, the ellipsoid distribution does not seem to capture the direction properly. Furthermore, it appears to completely rotate the principal axes in the air between the workpiece and an inductor by 90 • . It also drastically changes the direction of the axes of the ellipsoid on the surface of the workpiece. In contrast, the metric symmetrized by the T 2 procedure predicts a more uniform behavior. In fact, the ellipsoids point uniformly in  the same direction. Therefore, the T 2 procedure shows to be better suited to enforce symmetric and positive definiteness of the metric. Fig. 7 shows the resulting meshes after computing the metric tensors from the original isotropic mesh (Figs. 5 and 6) as well as their respective error estimator distribution over the domain. As predicted by the ellipsoid analysis, the highest element density is located over the inductor region, the workpiece close to the inductor, and the air in between, for both metrics. However, the mesh created from the T 1 procedure shows a change of orientation in the elements between the inductor and the workpiece, as predicted by the ellipsoids and a larger error compared to T 2 . Therefore, it is evident that procedure T 2 allows to generate a better suited mesh distribution, which leads to a higher decrease of the error.
In conclusion, following the results shown in Fig. 6 and 7, the T 2 procedure seems to be the most suitable to be used. In fact, the positive-definite symmetric matrix found here maintains more accurately the direction of the anisotropy, thus it will be used from here onward to generate the adaptative mesh.

B. Anisotropic Behavior
Once the appropriated symmetrization procedure has been chosen, a second test is performed with both metrics. The anisotropy behavior is then studied to determine which metric captures the physical phenomenon more accurately. The metric is symmetrized using (36) and as in the previous test, the ellipsoid is used to understand how each metric captures anisotropy. For this purpose, the magnetic field ( ⃗ H ) is used to calculate both metrics. Then, the characteristic ellipsoids of M ⃗ H and M ⃗ H are shown in Fig. 8. Fig. 8 shows a remarkable difference in their anisotropy. Although both metrics point uniformly in the same direction, the shape of the ellipsoid varies, indicating a higher anisotropy in the M ⃗ H metric. To better see both metric performances, the mesh is computed in both cases. The resulting meshes are displayed in Fig. 9(a) for M ⃗ H and Fig. 9(b) for M ⃗ H . Following the results of these resulting meshes, it is evident that the metric calculated from a vector field M ⃗ H follows the sharp anisotropy of the field. Whereas M ⃗ H shows also an anisotropic behavior but to a lesser extent.  The previous case is a quasi-2-D case of a unidirectional field. To probe both metrics in a more complex 3-D scenario, an axisymmetric case of IH is tested (see Fig. 10).
A similar comparison to the previous case is made in a 3-D model having a cylindrical field distribution. By comparing the ellipsoid representing each metric, the anisotropic behavior contained in the metric before meshing is clearly seen. These results are shown in Figs. 11 and 12   isotropic behavior. To see this more clearly, both metrics are used to drive the remeshing procedure, resulting in the meshes shown in Fig. 13. Fig. 13(a) and (b) shows the meshes generated with M ⃗ H and M ⃗ H , respectively. In the same way, as provided for the results with the representative ellipsoids, M ⃗ H has a higher anisotropic behavior. Although M ⃗ H is able to follow where the mesh needs to be refined or de-refined, its anisotropy is somewhat small.
Results presented in this section lead to the conclusion that the symmetrization procedure from the conjugate product (36) allows not only to symmetrize and positively define the tensor, but also does not distort the principal direction  of anisotropy. Moreover, although metrics computed as the Hessian of a scalar field have been more extensively studied in the literature, the metric computed from the Jacobian of a vector field seems to better capture the anisotropy of physical phenomena. In this sense, the M ⃗ H procedure will be retained and applied to several cases in Section V and results compared with the ones obtained with M ⃗ H .

V. APPLICATION TO EPMS
To test the capabilities of both metrics defined in Section III, several cases based on electromagnetically coupled processes are presented in this section. To begin, the intrinsic anisotropy of the electromagnetic phenomena is highlighted. Hence, the advantage of anisotropic mesh adaptation. Along the same lines, a comparison is made between an isotropic and an anisotropic approach in a benchmark case. Then, the metric is tested to prove its capacity to deal with the electromagnetic skin effect automatically. Finally, three industrial cases are presented: an IH application for heat treatment in a wheelbearing piece, a magnetic pulse-forming (MPF) application in a circular sheet, and an EMS application with a more complex geometry. The magnetic field ( ⃗ H ) is used to construct and compare the proposed metrics for almost all cases. Only the skin effect test is performed with the electrical field ( ⃗ E) with the purpose to show that either field can be used. These   implementations are carried out within the electromagnetic module of Forge. 2

A. Intrinsic Anisotropy
The electromagnetic problem solved in full immersion finite-element analysis does not require uniform accuracy throughout the mesh; a uniform mesh wastes computational resources, which is often the case in an isotropic configuration. It is thus more efficient to use a procedure for enabling adaptive anisotropic remeshing. Regions where fields change in a more abrupt manner are refined, and regions with small gradients can be discretized with larger elements. The anisotropic structure plays a key role in tracking these fields; the most obvious case is the one where the fields are fully single-directional, and we can allow the elements to be larger  in the direction orthogonal to the field, which results in decreasing the number of elements in the mesh while retaining good accuracy. We can take the example of Fig. 14, where a cylindrical inductor and a part are shown, and the magnetic field lines are also shown to clarify the phenomenon. Thanks to the magnetic field lines, it is possible to observe common electromagnetic features that create anisotropy, such as skin and proximity effects, the former pushing the field toward the surface of the cylindrical part and the latter making the field larger in the vicinity of the inductor.
Appropriate mesh adaptation will allow refinement in these regions and the growth of elements in other regions. Note that if we opt for an isotropic structure, the number of elements will grow exponentially as the length of the cylinders grows. To check this, an isotropic and an anisotropic mesh have been performed to contain the same error in both meshes. The M ⃗ H metric has been used to capture the anisotropic behavior of the magnetic field. Fig. 15 shows both isotropic and anisotropic mesh adaptations. For the same error, the isotropic mesh has 405 000 • of freedom, while   the anisotropic mesh needs 140 000. It is worth mentioning that the anisotropic mesh obtained with M ⃗ H is successful in capturing the electromagnetic phenomenon mentioned above.  Both meshes have been tested to compare the convergence of the finite-element solution. The conjugate gradient method is used to solve the linear system of the electromagnetic problem using the Jacobi preconditioner. Table III shows the CPU time and iteration needed to solve the linear system. Therefore, we can conclude that in this case, the anisotropic structure helps to reduce the computational resources by almost three times and makes the computation 40% faster.

B. Isotropic Versus Anisotropic Performance Comparison (Model Size/Error Distribution/Etc.)
To test the metric tensor developed in this work and compare it with an isotropic mesh adaptation procedure, a benchmark case has been selected. This benchmark case belongs to the testing electromagnetic analysis methods (TEAM) problems, which represents an open international working group aiming to compare electromagnetic analysis computer codes. Particularly, case number 7 [34]-more specifically focused on IH processes-has been selected. This problem consists of a thick aluminum plate with a hole located off-center and unsymmetrically in a nonuniform magnetic  field, which is used to construct the metrics. The field is produced by a sinusoidal current. It is shown in Fig. 16. Fig. 17 shows the mesh and error estimator of the comparison between isotropic and anisotropic mesh adaptation as well as the initial mesh used. The parameters to delimit the mesh size are defined the same in both cases (h min and h max ). As the space discretization used is a fully immersed finiteelement approach, it means that the whole domain is meshed. Different minimum and maximum mesh sizes are defined for each one of them (workpiece, inductor, and air). A frequency of 50 Hz and a current of 2742 are used. It can be seen in Fig. 17 and in Table IV that the isotropic approach, although it succeeds in following the physical phenomenon, the amount of degrees of freedom is almost twice as much as used for the anisotropic approach. It must be said that the anisotropic approach adequately follows the physical phenomenon and manages to distribute the error in the piece, even including several elements on the thickness of the plate. The error is well distributed in each of the approaches, however, the metric M ⃗ H achieves a better compromise between the error and degree of freedom.

C. Skin Effect
It is assumed that the generality of the metric tensors developed in this work can be adapted to any physical phenomenon. As explained above, the skin-depth phenomenon is one of those responsible for the anisotropy of the electromagnetic problem. A simple case to demonstrate this is presented in Fig. 10. The main idea is to vary the frequency parameter to obtain a different skin depth; using relation (11), frequencies of 350, 1000, and 550 Hz were used to obtain a skin depth of 12, 7, and 3 mm, respectively. A priori knowledge of the skin depth gives us an idea of whether the metric can capture this phenomenon and whether it manages to generate a new mesh to address it. In this case, only the metric calculated from the Jacobian is used using the electric   Fig. 18. It is worth mentioning that the radius of the workpiece is 15 mm. Fig. 18 shows the generated meshes for each skin depth parameter. It is clearly seen that the metric can adequately follow this phenomenon. Allowing a higher density of elements in the interior in the skin depth area and de-refining where it is no longer needed. Additionally, Fig. 19 shows the evolution of the mesh size on the horizontal axis with respect to the electric field and the error estimator. As these results reveal, it seems to be a good agreement between the fields. In particular, as the gradient of the electric field increases, the error estimator increases, and the metric M ⃗ E is naturally able to capture these changes and allows the mesh size to be decreased to try to capture the field more faithfully.

D. IH Application
The metric is tested in the case of IH by surface heat treatment of a wheel bearing. For this purpose, two inductors are placed in the area where the heat treatment is to be performed, and a current concentrator is placed on each inductor to direct the magnetic field toward the part. By taking advantage of the wheel bearing symmetry only a part of the geometry is simulated. The magnetic field ( ⃗ H ) is used to construct. At first glance, it is expected that the metric tensor tries to adapt the mesh to the regions between the part and the inductors where the field acts strongly, as well as capturing the special behavior added for the concentrators. The wheel bearing is shown in Fig. 20 and the results of a comparison of experiments versus simulation of induction hardening for the prediction of the microstructure is shown in Fig. 21, where the area of highest influence of the treatment is clearly seen [35]. A frequency of 13 020 Hz is used for the heat treatment. The results for this case are presented in Fig. 22. At first sight, it can be seen that both metrics are capable to follow the predicted electromagnetic behavior. Being able to properly refine the mesh on the bearing surface and making it coarser on the rest of the part, along with capturing the functionality of the concentrators. In addition to reducing the error on the adapted meshes. Although the metric M ⃗ H seems to  achieve a mesh with lower error than M ⃗ H , the number of degrees of freedom required in the latter case is lower as shown in Table V.

E. MPF Application
The MPF technology or electromagnetic forming (EMF) consists in deforming metallic components through the application of an intense electromagnetic pulse. This process enters the category of high-speed forming processes because of the range of strain rates that are usually attained, ranging from 10 3 to 10 4 s −1 . A typical coil and workpiece setup of this process is shown in Fig. 23. The objective of this case is to check if the metric tensor can capture the magnetic field induced from the circular coil at a high frequency, and how it manages to adapt the mesh. A discharge pulsed signal with a natural frequency of 162 kHz is used in this case and the magnetic field ( ⃗ H ) is also used to construct the metrics. Fig. 24 shows the comparison of the meshes and error distribution for both metrics and the initial mesh. As in the previous case, both metrics can accurately capture the physical phenomenon. It is again observed that the metric M ⃗ H achieves a smaller final error distribution than M ⃗ H . However, the difference between the degrees of freedom, shown in Table VI, is quite high in both of them, leading the metric M ⃗ H to require less degrees of freedom to achieve a mesh with an acceptable diminution of the error.

F. EMS Application
So far, metrics have been used to govern the anisotropic meshing procedure using simple geometries. To test its performance with more complex geometries, the EMS process has been chosen. This process is widely used in the materialsforming industry. It consists of using the time-varying electromagnetic field to control the fluid flow between the liquid steel and the stirrer without any physical contact. This process is used to disrupt the molten metal's fluid flow by means of the Lorentz force provided by a linear induction motor allowing a more homogeneous solidification, avoiding premature growth of dendrites during the casting process and resulting in a better quality of the final ingot. The geometry used to model this process was a steel bar and six inductors supported by a core, as shown in Fig. 25.
A low frequency characterizes this process, here 2 Hz has been used and the magnetic field ( ⃗ H ) is again used to construct the metrics. Fig. 26 shows the results for the EMS case. This case is an appropriate test to measure the performance of the metrics for dealing with complex geometries. The meshes and error distribution are also displayed as in the previous cases to complete the comparison between both proposed metrics. It is observed that the anisotropy is again more evident in the mesh generated using the metric M ⃗ H precisely on the workpiece. This leads to a mesh with less degrees of freedom, as shown in Table VII. It is worth mentioning that this time the metric M ⃗ H with smaller number of elements than M ⃗ H is able to achieve a better distribution of the error on the entire domain.

G. Magnetic and Electric Fields Duality
In the electromagnetic phenomena, we can find the duality between the magnetic and the electric fields. Therefore, leading to write the Maxwell equations in a simplified manner Since the metric proposed in this work allows working with any field. A comparison of the resulting meshes of ⃗ H and ⃗ E shows the duality between these fields in Fig. 27. In fact, constructing the metric of either of these fields generates an adapted mesh that follows the other field. Fig. 27(a) shows the mesh constructed from the electric field ⃗ E and the magnitude of the magnetic field and its vectors fields are shown. While Fig. 27(b) shows the mesh constructed from the magnetic field ⃗ H and the magnitude of the electric field and its vector fields are shown. It is seen in booths cases that the mesh is correctly adapted from their dual field.

VI. CONCLUSION
Several strategies for constructing a metric tensor to deal with electromagnetic phenomena have been presented. On the one hand, we showed two ways of symmetrizing the tensor field. The conjugate product is the most appropriate method to ensure symmetry and positive definiteness of the tensor, which also allows us to keep the direction of the anisotropy contained in the metric. Second, we showed that the metric field can be obtained as the Hessian of a scalar field or as the Jacobian of a vector field. Both procedures use a recovery-based method to apply the gradient operators. They have been tested in different configurations, from comparison with an isotropic mesh procedure to checking how they manage to capture the skin effect phenomenon, being sufficiently accurate to guarantee a level of refinement where is necessary. Industrial applications in induction heating, MPF, and EMS processes have then been presented, being able to cope with the physics involved in each of these processes, allowing good distribution and error reduction. This work has shown that the metric approximation as the Jacobian of a vector field could be a better option to carry out the mesh adaptation procedure, allowing to better capture the anisotropy of the fields and keeping a good balance between the error reduction and the number of degrees of freedom. It also provides an additional way to deal with the average anisotropy of a vector field. A more classical way would be to deal with the anisotropy of each component by computing a metric for each component and then carrying out metric intersection procedures. However, these procedures are more expensive and are also dependent on the order in which the intersections are carried out. This work has focused only on the electromagnetic model, leaving for the future the perspective of implementing it in a fully coupled multiphysical application.