thggnum.html

有限元学习研究用源代码(老外的),供科研人员参考

<h1><a NAME="numrequire">Requirements of the Numerical Method</a></h1>


In this section we want to consider the connection between the used
discretization scheme and the grid quality criteria we have to use.

<h2><a NAME="FEM-FV">Finite Elements versus Finite Volume</a></h2>


  <p>Let's consider the simplest case of the diffusion equation --- the
Laplace equation --- which describes the stationary state of a single
species with constant diffusion coefficient:
\Delta u = 0

Let's compare now the discretization of this standard example
operator for the two most popular discretization schemes --- the
Finite Element Method (FEM) and the Finite Volume (FV) discretization.
We consider only one question - the M-matrix property for the
resulting matrix.

  <p> The M-matrix property depends on the sign of the so-called
connection value between two neighbour nodes i and j. For FEM, this
connection value is defined by:
 \Gamma_{ij} = - \int \nabla \psi_i \nabla \psi_j dV

Here \psi_i is the basic function for the i-th node of the
discretization.
  <p>In the FV method, for every grid node a box around the node has to
be defined. For the most usual FV variant the Voronoi boxes are used.
We obtain the following connection value:
 \Gamma_{ij} = A_{ij} / l_{ij}

Here A_{ij} is the area of the common boundary of the Voronoi
boxes of the nodes i and j, and l_{ij} the distance between these
nodes. In above cases, the (weak) M-matrix property will be fulfilled
if

 \Gamma_{ij} >= 0

for every pair of neighbour nodes i and j.

  <p>In 2D, we have a satisfactory situation. The above coefficients
coincide, and we have the following well-known
 For a given set of points, there is only one grid --- the
Delaunay grid --- so that for every inner connection we have
\Gamma_{ij} >= 0 as for FEM, as for FV.  

  <p>It is well known that the generalization of this theorem to 3D
fails. The Delaunay grid is not the optimal grid for FEM.  Lo <a HREF="thbib.html#Lo1991"> [Lo1991] </a> thinks that this is a fundamental weakness
of 3D Delaunay triangulation.  The typical situation is the so-called
sliver --- a badly distorted tetrahedron with well-proportioned faces
but with arbitrary small volume. Sliver can account for as much as 10
% of the total number of tetrahedra generation. They lead to
arbitrary big terms with incorrect sign of the related connection
value \Gamma_{ij}.

  <p>Sometimes sliver can be removed by a local grid
transformation. But, in general, this does not lead to the M-matrix
property for the grid. Letniowski <a HREF="thbib.html#Letniowski1992"> [Letniowski1992] </a>
has given a simple example of a point set so that every grid has an
inner connection with \Gamma_{ij} &lt; 0.

  <p>Our opinion is slightly different from the point of view of Lo.
If his conclusion is "How bad for the Delaunay grid", we say "How bad
for the FEM method". The reason is that the generalization fails only
for FEM, not for FV:
 For the FV method, for a given set of points, there is only
one grid --- the Delaunay grid --- so that for every inner connection
we have \Gamma_{ij} >= 0.  

  <p>The proof is a trivial consequence of the formula we have given before
for the FV method.

  <p>Remark that for the time-derivative part of the diffusion equation\partial_t u very often "mass matrix lumping" will be used because
of the numerical problems of the pure FEM formula. This modification
also may be considered as going from FEM to some FV method. There
seems to be no reason not to do the same for the Laplace term.

<h2><a NAME="nonconstdiff">Non-Constant Diffusion Coefficient</a></h2>


Consider now shortly how we have to handle a non-constant diffusion
coefficient.

  <p>One, very popular point of view is to consider the diffusion
coefficient as constant over each element. But this strategy can
easily lead to an incorrect sign for a connection coefficient.  Remark
that the connection coefficient \Gamma_{ij} may be described as
the sum over all elements containing these two nodes. The previous
theorem guarantees that this sum is greater zero for the Delaunay
grid, but not that every part is greater zero. If these parts will be
weighted in another way by the different diffusion coefficients of the
elements, we can obtain resulting coefficients with incorrect sign.

  <p>If we consider the diffusion coefficient as constant over the
"domain of influence" of the edges (consisting of the two cones over
the common boundary of the Voronoi cells), we can avoid this effect.
In this case, we have to compute the diffusion coefficient over this
domain and to multiply it with the coefficient \Gamma_{ij}.

<h2><a NAME="boundconsid">Boundary Consideration</a></h2>


For boundary edges the previous theorems do not lead to a correct sign
of the connection coefficient \Gamma_{ij}. Here we have to modify
the point set itself to obtain a grid which leads to an M-matrix
property.

  <p>The standard technique in this case is to include new boundary
points. Usually this technique allows to avoid connection values with
incorrect sign. In the case of anisotropic grids, this may be not
the ideal solution, and other techniques have to be used to obtain
this result.

  <p>If a point set fulfills the M-matrix condition also for the
internal boundaries, also another problem of Delaunay grids will be
solved. In general, the Delaunay grid may not preserve internal
boundaries. There may be Delaunay edges between internal nodes of
different regions. But we have the following
If the grid consisting of different sub-grids for different
regions leads to an M-matrix property for the FV method separately for
each sub-grid, it will be the Delaunay grid of this point set.

Especially, in this case the Delaunay grid preserves the inner
boundaries.

<h2><a NAME="ggnumconclusion">Conclusion</a></h2>


We see that we can the M-matrix property by using the following
strategy:

 
 <ul>
 <li>Using Delaunay grids.
 
 <li>Using Finite Volume instead of FEM.
 
 <li>Using averages over edges instead of elements for diffusion
coefficients.
 
 <li>Using point sets good enough near the boundaries. Especially the
Delaunay grid must preserve the boundaries.
 
 </ul>


One question is if the M-matrix property is really so important.  An
M-matrix is desirable for iterative sparse matrix methods. But it is
also possible to solve problems without this property.  The main
reason for us is that with an M-matrix the approximation satisfies a
discrete maximum principle. This leads to nonnegative concentrations
if the initial conditions and boundary conditions are correct. Remark,
that this is very important especially for diffusion equations,
because negative concentrations are physically nonsense. For other
applications of the Laplace operator the importance may be not so big.
So, for other applications higher accuracy combined with some small
negative values may be the better solution.