{'fraction_non_alphanumeric': 0.03574729337206232, 'fraction_numerical': 0.019804594665962502, 'mean_word_length': 4.329287598944591, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Anisotropic mesh adaptation for 3D anisotropic diffusion problems with application to fractured reservoir simulation

Anisotropic mesh adaptation is studied for linear finite element solution of 3D anisotropic diffusion problems. The M-uniform mesh approach is used, where an anisotropic adaptive mesh is generated as a uniform one in the metric specified by a tensor. In addition to mesh adaptation, preservation of the maximum principle is also studied. Some new sufficient conditions for maximum principle preservation are developed, and a mesh quality measure is defined to server as a good indicator. Four different metric tensors are investigated: one is the identity matrix, one focuses on minimizing an error bound, another one on preservation of the maximum principle, while the fourth combines both. Numerical examples show that these metric tensors serve their purposes. Particularly, the fourth leads to meshes that improve the satisfaction of the maximum principle by the finite element solution while concentrating elements in regions where the error is large. Application of the anisotropic mesh adaptation to fractured reservoir simulation in petroleum engineering is also investigated, where unphysical solutions can occur and mesh adaptation can help improving the satisfaction of the maximum principle.


Introduction
We are concerned with the linear finite element solution of the three dimensional boundary value problem (BVP) of the diffusion equation in Ω (1) subject to the Dirichlet boundary condition u = g, on ∂Ω (2) where Ω is a bounded polyhedral domain, f and g are given functions, and D is the diffusion matrix. We assume that D = D(x) is a symmetric and uniformly positive definite matrix-valued function on Ω. It includes both isotropic and anisotropic diffusion as special examples. In the former case, D takes the form D = α(x)I, where I is the 3 × 3 identity matrix and α = α(x) is a scalar function. In the latter case, on the other hand, D has not-all-equal eigenvalues at least on a portion of Ω. Anisotropic diffusion problems arise from various branches of science and engineering, including plasma physics [21,22,23,51,57], petroleum engineering [1,2,15,50], and image processing [10,11,34,54,62]. When a conventional numerical method is used to solve the problems, spurious oscillations may occur in computed solutions. Numerous research has been done for two dimensional (2D) problems; among other works, we mention a few here, [8,12,18,23,32,33,38,40,41,42,43,44,55,61,63,65,66]. A common approach is to design a proper discretization method and/or a proper mesh so that the numerical solution satisfies the maximum principle (MP). Recently, an anisotropic non-obtuse angle condition was developed in [31,40,41,46] for the linear finite element solution of both time independent and time dependent anisotropic diffusion problems to satisfy MP.
On the other hand, much less work has been done for 3D anisotropic diffusion problems. Although MP preservation has been studied in general dimensions e.g. in [8,12,23,33,38,37,40,41], most of them either consider isotropic diffusion or present numerical examples only in 1D and 2D. For example, only isotropic diffusion is considered in [12,37]. It is shown in [39] that the 3D Delaunay triangulation does not generally produce a mesh with which the numerical solution satisfies MP. Mesh conditions are studied in [8] for a reaction-isotropic-diffusion problem for general dimensions and numerical examples in 1D and 2D are presented. The difficulty of MP satisfaction for 3D problems is remarked in both [8] and [39].
The objective of this paper is to study the linear finite element solution of 3D anisotropic diffusion problems. The focus will be on MP preservation and mesh adaptation. Four different metric tensors used in anisotropic mesh generation will be considered. This study is a 3D extension of the work [40]. Moreover, new sufficient conditions will be developed for the linear finite element approximation to satisfies MP, and a mesh quality measure will be defined to provide a useful indication for MP satisfaction. The mesh quality measure is developed along the approach of [28]. But the interested reader is also referred to [35,36] and references therein for different mesh optimization methods and quality metrics. Furthermore, the application to fractured reservoir simulation in petroleum engineering will also be investigated, where unphysical solutions can occur and mesh adaptation can help improving MP satisfaction.
An outline of the paper is given as follows. The linear finite element formulation for BVP (1) and (2) is given in Section 2. MP preservation and some sufficient conditions will be discussed. Section 3 contains the discussion of anisotropic mesh adaptation and four metric tensors. Numerical examples are given in Section 4, followed by the investigation of application of anisotropic mesh adaptation to fractured reservoir simulation in Section 5. Conclusions are drawn in Section 6.

Finite Element Formulation
In this section we briefly describe the piecewise linear finite element discretization for BVP (1) and (2) and state a few properties of the discretization.
Let U g = {v ∈ H 1 (Ω) | v| ∂Ω = g}. The weak formulation of BVP (1) and (2) is to find u ∈ U g such that For the finite element discretization, we assume that a tetrahedral mesh T h is given on Ω. Let g h be the piecewise linear interpolation of g on the boundary ∂Ω and U h g h be the piecewise linear finite element space on T h with the boundary data g h . Then, the finite element formulation for BVP (1) and (2) This discretization is standard and it is expected that u h converges to u at a rate of second order in the L 2 norm and first order in H 1 norm. We take the reference elementK to be a unitary equilateral tetrahedron and denote an element in the mesh T h by K. Moreover, denote byq i (i = 1, ..., 4) the unit inward normal to the face facing the i th vertex ofK. Let F K be the affine mapping fromK to K and F K the Jacobian matrix of F K .
Theorem 2.1 (Li and Huang [40]). If the mesh satisfieŝ where D K is the average of D over K, then the linear finite element scheme (3) for solving BVP (1) and (2) satisfies the maximum principle, Preserving MP is crucial to avoiding artificial oscillations in the computed solution. It is thus interesting to know what meshes satisfy the anisotropic non-obtuse angle condition (4). Obviously, a sufficient condition is where c K is a positive scalar constant on K and I is the 3 × 3 identity matrix. A weaker condition is stated in the following theorem, where a general case for the d space dimension (d ≥ 2) is considered.
Notice that the bound in (7) has the value Moreover, (5) is sufficient for (7) because (5) implies that all of the eigenvalues of (F K ) −1 D K (F K ) −T are equal to c K and thus the left-hand side of (7) is equal to one which is less than the right-hand side. This means that (7) is weaker than (5).
Unfortunately, (7) is still stronger than (4). Consider the triangular and tetrahedral elements in Fig. 1 and the case with D = I. It is easy to see that the elements are non-obtuse and satisfy (4). A direct calculation shows that 1.732, for 2D elements in Fig. 1(a) 2.151, for 3D elements in Fig. 1(b) which violates (7). Nevertheless, as will be seen in the next section, the conditions (5) and (7) provide useful guidelines for generating meshes that improve the MP satisfaction of the scheme.

Anisotropic mesh adaptation
In this section we study anisotropic mesh adaptation for the anisotropic diffusion problem (1) and (2). We use the so-called M-uniform mesh approach [29,30] where an adaptive mesh is viewed as a uniform one in the metric specified by a tensor M = M(x) which is assumed to be symmetric and uniformly positive definite on Ω. For the moment, we assume that M has been given. Its choice will be discussed later. It is shown in [30] that an M-uniform tetrahedral mesh T h satisfies where Condition (10) is called as the equidistribution condition which requires all of the elements to have the same size in the metric M K and therefore determines the size of K from det(M K ) 1 2 . On the other hand, (11) is called the alignment condition which requires K, measured in the metric M K , to be similar to the reference elementK that is measured in the Euclidean metric and thus controls the shape and orientation of K. It can also be shown that the principal axes of the circumscribed ellipsoid of K are required to be parallel to the eigenvectors of M K while their lengths are reciprocally proportional to the square roots of the respective eigenvalues [30].
M-uniform or nearly M-uniform meshes can be generated using various strategies; see [30, Section 4] for more detailed discussion. We mention just a few here, the Delaunay triangulation method [5,6,9,53], the advancing front method [20,48], the bubble mesh method [64], the combination of combining refinement, local modification, and local node movement [3,7,17,24,45], and the variational method [27]. Particularly, we mention two C++ codes which can directly take the user supplied metric tensor. One is BAMG (Bidimensional Anisotropic Mesh Generator) for 2D meshes developed by Hecht [25] based on the Delaunay triangulation and local node movement. The other is MMG3D (Anisotropic Tetrahedral Remesher/Moving Mesh Generation) developed by Dobrzynski and Frey [16] based on refinement and local node movement. The latter is used in our computation.
A key component of the M-uniform mesh approach is to define the metric tensor. We consider four choices here. The first one is This is the simplest choice, and the resulting meshes are uniform or nearly uniform. The second choice is based on linear interpolation error. It is known that the error for the linear finite element solution to (1) and (2) is bounded by the error in the piecewise linear interpolation of the exact solution. Thus, it is reasonable to define the metric tensor based on linear interpolation error. A metric tensor based on minimization of the H 1 semi-norm of linear interpolation error is given [29] by where u h is the finite element solution, H K (u h ) is a recovered Hessian of u h over K, |H K (u h )| is the eigen-decomposition of H K (u h ) with the eigenvalues being replaced by their absolute values, and α h is defined implicitly through This equation determines α h uniquely and can be solved, for instance, by the bisection method. Moreover, the choice of α h in this way concentrates roughly 50% of the mesh points in the region where the solution changes significantly. The third and fourth choices are related to the diffusion matrix. Recall that (11) requires that all of the eigenvalues of the matrix ( for some constantc K on K. Comparing this with (5) suggests that we choose the metric tensor as where θ K is an arbitrary positive piecewise constant function and D K is the average of D over K. From (16), it is not difficult to see that any M-uniform mesh associated with this metric tensor satisfies (5) and therefore, the finite element solution u h satisfies MP. Then, the third choice is which corresponds to θ K = 1. For the fourth choice, we take advantage of the arbitrariness of θ K in (17) and choose it to minimize a bound of H 1 semi-norm of linear interpolation error. This way we can combine the desire of preserving MP with mesh adaptation. The metric tensor reads as (cf. [40]) where The M-uniform meshes associated with M DM P and M DM P +adap satisfy the alignment condition (11), and the mesh elements are aligned along the principal diffusion direction defined as the direction of the eigenvector corresponding to the largest eigenvalue of the diffusion matrix D.
It should be pointed out that in practical computation, it is difficult, if not impossible, to generate a perfect M-uniform mesh that satisfy the conditions (10) and (11). Thus, it makes sense to measure how far a given mesh is from satisfying these conditions. From (10) and (11), we can define the equidistribution and alignment measures as It is not difficult to show that the maximum norm of both functions, Q eq L ∞ and Q ali L ∞ , has the range [1, ∞). Moreover, Q eq L ∞ = 1 and Q ali L ∞ = 1 imply a perfect M-uniform mesh. The larger Q eq L ∞ and Q ali L ∞ are, the farther the mesh is away from being M-uniform. It is worth emphasizing that only Q ali is related to MP satisfaction. The quantity Q eq indicates how evenly the error is distributed among the elements. It is interesting to point out that the definition (21) is slightly different from a more common definition [30] where the trace of ( These two definitions are equivalent since the trace and 2-norm of a positive definite matrix are equivalent. The main motivation for which we use (21) is that the condition (7) can now be rewritten as As mentioned before, this condition is hard to satisfy. Indeed, for meshes shown in Fig. 1, we have Q ali L ∞ = 1.732 (2D) and 2.151 (3D). Although they violate (23), the meshes can still be considered to be close to satisfying the alignment condition since Q ali L ∞ is relatively small. Moreover, (23) shows a connection between the alignment condition and the preservation of MP.

Numerical examples
In this section we present two three-dimensional examples to demonstrate the performance of the anisotropic mesh adaptation strategy described in the previous section with the four metric tensors. An iterative procedure for solving PDEs using anisotropic adaptive meshes is shown in Fig. 2. The PDE is first solved on the current mesh using piecewise linear finite elements. Then, the metric tensor M is computed, followed by the generation of a new mesh for the metric tensor using MMG3D. This procedure is repeated a few times untill the quliaty measures (20) and (21) cannot not be improved further. In our computations, the procedure is repeated five times and numerical experiment shows that there is no significant improvement in the results when more iterations are used. The computations are performed in the framework of an open source finite element software, FreeFem++ developed by Hecht [26], with the linear solver being chosen as the conjugate gradient method with tolerance 10 −15 . Moreover, the Hessian of the finite element solution u h used in the computation of the metric tensors is computed using the "mshmet" library embedded in FreeFem++. The "mshmet" library first approximates the gradient (first derivatives) using the least squares fitting, then uses the approximated gradient in the least squares fitting for the Hessian matrix (second derivatives).

Given a mesh
Solve PDE

Recover solution derivatives and compute M
Generate new mesh according to M The meshes associated with the four metric tensors will be denoted with the same notation for the metric tensors. For instance, a mesh associated with M adap will be called an M adap mesh. We consider the diffusion matrix in the form where k 1 is the dominant eigenvalue, φ is the angle between the principal diffusion direction and the positive z-axis, and θ is the angle between the projection of the principal diffusion vector on the xy-plane and the positive x-axis.
Example 4.1. For the first example, we would like to compare the L 2 -norm of solution errors obtained from different meshes and check the violation/satisfaction of MP. We consider problem (1) in the unit cube Ω = (0, 1) 3 . D is chosen in the form of (24) with φ = −π/4, θ = 5π/6, k 1 = 100, k 2 = 10, and k 3 = 1. The source function f and the boundary function g are chosen such that the exact solution is given by The continuous problem satisfies MP and the solution is in the interval (0, 1 + e]. The exact solution in the domain Ω and on the cross-sections x = 0.5 and z = 0.5 is shown in Fig. 3. Fig. 4 shows the four types of meshes with a view of the inner mesh by cutting a corner of the domain. Fig. 5 shows the meshes on the cross-section x = 0.5. The elements in the M DM P mesh are aligned well along the primary diffusion direction θ = 5π/6, while the elements in the M id mesh are aligned along the direction of θ = π/4. The elements in the M adap mesh concentrate around the central region where the solution changes rapidly. The elements in the M DM P +adap mesh not only are aligned along the primary diffusion direction but also concentrate around the central region, which is a result of combining MP preservation and adaptivity. Table 1 shows the L 2 -norm of the solution error obtained on M id , M adap , M DM P , and M DM P +adap meshes. The minimal value in the computed solution, which is denoted by u min and indicates the undershoots and violation of MP, is also reported in the table. The convergence history of the L 2 -norm of the error is shown in Fig. 6. It is clear that the error converges at a rate of second order as the mesh is refined. Moreover, the numerical solution obtained using the M adap mesh has the smallest error, which is consistent with the fact that M adap is formulated based on the minimization of the interpolation error (and the finite element error). The solution obtained using the M DM P mesh has the largest error due to the fact that the mesh is designed to satisfy MP, which is generally in conflict with error reduction. The error of the solution obtained using M DM P +adap is between those using M adap and M DM P , which is expected from the construction of M DM P +adap and shows a good balance between mesh adaptivity and MP preservation.
It can also be observed from Table 1 that the numerical solution obtained from M id and M adap violate MP even for very fine meshes. On the other hand, the numerical solutions obtained from M DM P and M DM P +adap violates MP for coarse meshes but satisfy MP when the mesh is sufficiently fine. This indicates that M DM P and M DM P +adap can help improving the satisfaction of MP for the finite element solution.  We first consider two different values of θ and then a diffusion matrix with a strong anisotropic feature. The continuous problem satisfies MP and the exact solution satisfies 0 ≤ u ≤ 4. Fig. 7 shows the numerical solution obtained using a fine M id mesh of 5,952,000 tetrahedra. Four different meshes are shown in Fig. 8. Table 2 shows the minimal value u min in the numerical solution obtained using different meshes. The corresponding mesh quality measures ||Q ali || L 2 , ||Q ali || L ∞ , and ||Q eq || L 2 are also reported in Table 2. Recall that only Q ali is related to MP satisfaction, while Q eq indicates how closely the mesh satisfies the equidistribution condition (10). The mesh quality measures for M id mesh are calculated using M = D −1 in order to see if it satisfies (23). The numerical solutions from all the meshes have the maximal value of u max = 4. It can be seen that the numerical solution obtained using the M adap mesh has undershoots (u min < 0) while the solutions obtained using other meshes satisfy MP. Interestingly, for this case MMG3D produces the same mesh from M DM P and M id . An explanation is that the initial mesh, the M id mesh, is already very close to an M DM P mesh (with ||Q ali || L 2 = 2.32). Since a half of the M id mesh elements are aligned with the primary diffusion direction and the anisotropy is not significant, the elements do not need to be stretched more. In fact, the M id mesh is closer to an M DM P mesh when k 1 = 50 with ||Q ali || L 2 = 1.72. When the anisotropy is more significant, the M id mesh is farther away from the M DM P mesh and the elements need to be stretched more along the primary diffusion direction. MMG3D will then adapt the mesh to be more different from the M id mesh, as will be shown in Case 3. Moreover, both M adap and M DM P +adap meshes have elements of worse shape (with large ||Q ali || L ∞ ) but the overall quality is acceptable (with relatively small ||Q ali || L 2 ). Case 2. θ = 3π/4. In this case, the elements of the M id mesh are no longer aligned along the primary diffusion direction. Meshes on the cross-section of z = y are shown in Fig. 9. Notice that the orientation of elements in the M DM P mesh (3π/4 direction) is very different from that in the M id mesh (π/4 direction). The results of u min and mesh quality measures are listed in Table 3.
One can see that none of the meshes satisfies (23), however, recall that (23) or (7) is only a sufficient condition for the MP satisfaction. Notice that the numerical solutions obtained from both the M DM P mesh and M DM P +adap mesh satisfy MP whereas others do not. Fig. 9 also shows that M DM P +adap mesh not only tempts to make elements to be aligned with the primary diffusion direction but also concentrates elements to minimize the solution error.
Case 3. θ = π/4, k 1 = 1000, k 2 = 1, k 3 = 1. In this case, the diffusion is much faster in the direction of θ = π/4 than in other directions. Fig. 10 shows the numerical solution using a fine M id mesh with 5,952,000 tetrahedra. A planar view of four different meshes at cross-section z = y is shown in Fig. 11. The results of u min and mesh quality measures are reported in Table 4.
Overall, all but the M id mesh are close to being M-uniform, with relatively small Q ali L 2 and Q eq L 2 . But they do not satisfy (7) so there is no guarantee that the solutions will satisfy MP.  Table 4

Application to fractured reservoir simulation in petroleum engineering
Numerical simulation plays an important role in petroleum engineering to predict production rate, optimize hydraulic fracturing design, and evaluate enhanced oil recovery processes. In those computations, the mesh has to be sufficiently refined around wellbore or fractures to accurately represent the flow effects thereon because the fractures are significantly smaller but with much higher permeability comparing to the reservoir matrix. It is challenging to consider the mesh inside a fracture due to its small scale. In recent years, shale gas reservoirs have become an attractive source for natural gas production largely due to the advancement of horizontal drilling and hydraulic fracturing techniques [60]. Shale reservoirs typically have extremely low permeability. For example, the representative permeability is 1.0 × 10 −4 millidarcies (mD) in the Barnett shale but is 5.0 × 10 4 mD in fractures [52]. The high permeability in fractures makes it possible to produce the gas from the shale while it makes the reservoir highly heterogeneous and anisotropic. The complexity of fracture network together with the complexity of shale pore structure makes shale-gas reservoir simulation a very challenging task. Some researchers focus on fluid flow models in nano-scale shale pores [13,47,49,58] and some other focus on different models of fracture network in the reservoir [14,19,52,56].
As our first exploration of the use of anisotropic mesh adaptation in fractured reservoir simulation, we consider the steady-state fluid flow that can be applied to effective permeability calculation [4,59]. For simplicity, we consider incompressible single-phase fluid flow and choose the single porosity dual permeability model. Darcy's Law is chosen to describe the flow in both matrix and fractures with different effective permeability in the corresponding regions. We also assume that permeability does not change with pressure. Under this setting, we can perform mesh adaptation in both the matrix We consider a reservoir with 3000 ft in x-direction, 1000 ft in y-direction, and 300 ft in z-direction, and denote the domain as Ω. Fig. 12 shows the sketch of the half-panel of the reservoir (0 ≤ y ≤ 500). A horizontal well is located in the center of the reservoir with length L w = 1400 ft along x-direction. The radius of the wellbore is 0.3 ft. The reservoir pressure is P r = 3800 psi and the pressure in the wellbore is P w = 1000 psi. There are three fractures with half-length L f = 400 ft and width W f = 0.01 ft at the location x = 800 ft, x = 1500 ft, and x = 2200 ft, respectively. The angle between the fractures and the positive x-axis is π/4. The height of the fractures is the same as the height of the reservoir. The subdomain formed by the fractures is denoted as Ω f . The permeability in the matrix is k mpar = 0.1 mD in the direction parallel to the fractures, is k mperp = 5 mD in the direction perpendicular to the the fractures on xy-plane, and is k mpar in the direction of z-axis. The permeability in the fractures is k pf = 10 4 mD along the fractures and is considered the same as that permeability in the matrix in other directions. For successful computation, a common technique is to scale up the width of the fracture while keeping the same conductivity, that is, K pf W f = 100 mD-ft for this example. In our computations, the width of each fracture is scaled up to W f = 10 ft and the permeability is reduced to K pf = 10 mD. Outside of the fractures, the primary diffusion direction in the matrix is perpendicular to the fractures.
The faces on the left side (x = 0), right side (x = 3000), and back side (y = 500) are denoted as Γ 1 . Γ 2 denotes the front side (y = 0), and Γ 3 denotes the bottom side (z = 0) and top side (z = 300). Γ 4 denotes the faces of the fractures on the front side Γ 2 . The horizontal well is in the center of the face Γ 2 along the x-direction. For simplicity, we assume that the pressure on the whole surface Γ 4 is the same as the pressure in the wellbore P w = 1000 psi. The pressure on Γ 1 is always the same as the reservoir pressure P r = 3800 psi. There is no flow through Γ 2 except in Γ 4 , and no flow through Γ 3 either. The mathematical formulation of the model is given by where the viscosity of the fluid and the porosity of the matrix are taken out from the equation from the assumption that they are independent of pressure. Gravitational effects are also ignored. The numerical solution with a fine mesh of 2, 804, 175 tetrahedra is shown in Fig. 13, where the mesh is manually refined around the fractures. The xy-planar view of four different meshes M id , M adap , M DM P , and M DM P +adap at cross-section z = 100 ft are shown in Fig. 14. In this example, the M id mesh is the initial mesh that has manual local refinement around the fractures. The numerical solutions obtained from M id and M adap meshes have pressure values larger than 3800 psi inside the domain, which violates MP and causes spurious back flow as shown in Fig. 15. The red shaded regions indicate unphysical pressure values. The numerical solutions obtained using M DM P and M DM P +adap meshes both satisfy MP, and the gradients of pressure are shown in Fig. 16. Different meshes at cross-sections along the fracture at x = 1500 ft and around the fractures at cross-section of z = 100 ft are shown in Figs. 17 and 18, respectively. M DM P and M DM P +adap meshes clearly are aligned along the principal diffusion direction both in fracture and in the matrix (where the primary diffusion direction is perpendicular to fractures). The alignment of mesh elements helps balance the anisotropic diffusion, and the numerical solutions obtained using both the M DM P and M DM P +adap meshes satisfy MP. M DM P +adap mesh also shows a degree of element concentration around the fractures. Table 5 lists the maximum pressure values for different meshes, which shows that improper meshes can lead to unphysical solutions in the fractured reservoir simulation.

Conclusions and comments
In the previous sections we have studied anisotropic mesh adaptation and maximum principle preservation for the finite element solution of three-dimensional anisotropic diffusion problems. We have defined a quality measure Q ali in (21) that provides a good guidance for MP satisfaction. And we have also developed some new sufficient conditions (cf. (5), (7), and (23)) that guarantee the linear finite element solution to satisfy MP. Comparison has been performed for four different metric tensors that are used to generate adaptive meshes with the existing software MMG3D. One of the metric tensors is the identity matrix which typically leads to uniform or almost uniform meshes. The other three are related to the diffusion matrix and/or the finite element solution. M adap (14) is based on minimizing H 1 semi-norm of linear interpolation error. M adap meshes concentrate elements near the regions where the Hessian of the solution is large, and gives the smallest error. M DM P (18) is taken as the inverse of the diffusion matrix and leads to meshes with elements aligned along the primary diffusion directions, which improves the MP satisfaction of the finite element solution but gives the largest error among all of the considered metric tensors. The last metric tensor, M DM P +adap (19) is proportional to the inverse of the diffusion matrix, with the coefficient of proportionality being The anisotropic mesh adaptation procedure is applied to fractured reservoir simulation in petroleum engineering. Numerical results show that M id mesh (with manual adaptation around fractures) can lead to unphysical solutions, M adap is capable of concentrating mesh elements around the fractures to reduce solution errors but may still violate MP, while both M DM P and M DM P +adap tend to align elements along the primary diffusion direction where the permeability is significantly larger than those in other directions and produce no artificial oscillations in the computed solution.
Overall, the numerical observations we made here for three dimensional problems are consistent with those for two dimensions reported in [40]. However, numerical experience suggests that it is much harder in 3D to make the mesh to satisfy or closely satisfy the alignment condition (11). This is especially true for the mesh adaptation situation (with the metric tensor depending on the solution) for which Q ali is relatively large for some elements although Q ali L 2 (which is an average of Q ali ) stays relatively small. How to generate better nearly M-uniform meshes for a given metric tensor in 3D has attracted interest from researchers (e.g. see [45,48]) and certainly deserves more investigations.   [25] F. Hecht. Bidimensional anisotropic mesh generator software (bamg). http://www.ann.jussieu.fr/hecht/ftp/bamg/bamg-v1.01.tar.gz, 2010.