thggdel.html

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

if the reverse test fails (f.e. if the target point lies on the
plane). A good neighbour is obviously allowed.

  <p>The idea of the algorithm is simple. We travel through the grid
by moving from the current element to any good neighbour until
we have found an element without any good neighbour.

  <p>If there is no allowed neighbour, we have really found the element
containing the point. If we have only one allowed neighbour, the
target lies in the neighbourhood of the side between them. It is easy
to establish a bound for the distance. If the target point lies in the
neighbourhood of an edge or a node, there may be two or three allowed
neighbours. In this case, the bound for the distance of the target
point from the element depends on the interior angles of the element.
But, independent of this distance, the algorithm has successfully
finished.

  <p>The question is if this algorithm always finishes. And really, it
is possible to construct some artificial grids where our algorithm can
lead to an infinite cycle.

  <p>A minor modification of the algorithm can avoid this. We count the
number of visits to each element. A neighbour which is allowed, but
not good, we include into the consideration with an initial penalty of
one.  If there are good neighbours, we go to one of the allowed
neighbours with the minimal count.

  <p>Now we can easily prove that the algorithm will not lead to an
infinite loop:

 The previously defined search algorithm ends in finite time
in an element which contains the target point in a small
neighbourhood. 

 Consider the line from an arbitrary point inside the
element. If this line intersects grid objects of codimension 2 and
higher, modify the start point slightly into general position.
Travelling along this line we obtain in exact arithmetics a sequence
of elements so that each is an allowed neighbour of the previous. If
we have an infinite loop, there will be elements with an unbounded
number of visits during the infinite loop. Consider the sequence from
such an element to the target element. Because the target element has
no visits (else the algorithm finishes because there are no good
neighbours), there must be an element with unbounded number of visits
with an allowed neighbour with bounded number of visits. This is not
possible.
 

What can we say about the algorithm without counting?

For a non-degenerated Delaunay grid, the search algorithm ends in
finite time even without the counting. 

If we go to a good neighbour in a Delaunay grid, the distance between
the centre of the circumscribing hyper-sphere becomes smaller.

  <p>In 2D, it is easy to find a (non-Delaunay) grid which allows an
infinite loop. Based on such an example, we can construct also a 3D
degenerated Delaunay grid with such an infinite loop. We simply have
to consider the grid consisting of the north pole of a big sphere and
the projections of this 2D grid from the north pole on this
sphere. Because all nodes of this grid are located on the same sphere,
every grid will be a degenerated Delaunay grid. Especially the grid we
obtain by connecting the north pole with the triangles of the given 2D
grid.

  <p>Thus, if we have Delaunay grids, there is de-facto no danger that the
search will lead to an infinite loop. Only for applications which require
very high security standards, the variant with counting seems necessary.
<h2><a NAME="ggdeltime">Time Requirements of the Algorithm </a></h2>


  <p> The Delaunay grid generation is the most time-consuming part in
our implementation. So, let's consider at first the time dependence
problem and our way to solve it. The efficiency of the algorithm
depends on different things:

 
 <dl>
 <dt>On the point set:
 
 <dd>In the worst case, the triangulation of n points contains
O(n^2) tetrahedrons. That's why, the general efficiency
cannot be better. But in our algorithm it is possible to avoid bad
situations (like a spider) in the step of point set generation so that
the number of tetrahedrons per grid point is bounded.

 
 <dt>On the point order:
 
 <dd>It is possible, that there is no spider in the resulting grid,
but they occur in intermediate steps. This intermediate creation and
removal of tetrahedra may require a lot of time. So, it is useful to
have this in mind implementing the point set generation. In our
algorithm, the order of point inclusion is the order of point creation
in the steps before. We try to avoid this effect in the organization of
the refinement step.

 
 <dt>On the search algorithm:
 
 <dd> The standard neighbourhood search algorithm requires O(n^{1 +
{1 \over n}}) for the search. Using tree search techniques it is possible
to reduce this time to O(n \log n). In our case, we can use
additional information from the octree data structure --- the
information about the nearest previously inserted point. This allows
search in linear time.
 
 </dl>


If our attempts to avoid the second effect are successful, this leads
in the result to an algorithm with linear behaviour. The empirical
results coincide with this dependence. Sometimes in numerical
experiments the algorithm looks even better then linear. The reason
seems to be that higher refinement leads in these cases to a more
regular grid.

<h2><a NAME="gdsBCorrect">Boundary Correction </a></h2>


  <p>The remaining problem connected with the Delaunay algorithm is
that this algorithm does not consider boundaries. So, the boundary may be
intersected by elements. The resulting grid is incorrect. We have the
following strategies to solve this problem:

 
 <ul>
 <li> The optimal strategy is to create a point set without such
effects. This is in good agreement with our considerations about
numerical quality. If we have a grid which fulfills the FV M-matrix
property for each domain, it will be the Delaunay grid of this point
set.
 
 <li> We can rebuild the grid using local transformations. In 3D this
will be very complicate. For 2D it can be proved that the swapping
procedure it sufficient to transfer every grid into a Delaunay
grid. In 3D the generalization of this algorithm can fail. The other
problem is that the resulting grid is not a Delaunay grid, so it's
numerical quality is not very good.
 
 <li> We can shift one point of the edge intersecting the boundary to
the intersection point. The problem is that around this point we
obtain in general a non-Delaunay or even degenerated grid. We have to
test this. Local transformations may be used to recover the Delaunay
property, but this can lead to a spider.
 
 <li> We can divide an edge intersecting the boundary at the
intersection point. Here we can also loose the Delaunay property or
obtain a spider.
 
 </ul>


We see that the boundary correction is a difficult problem. The best
way to solve this problem is the first. But there are situations where
the first method requires to a lot of additional nodes, or it has
failed because of a constellation we have not considered. In such
cases, the other methods have to be used as the "emergency exit" to
avoid a failure.

  <p>In the case of isotropic refinement, if the spider is not
dangerous, the last method seems to be universal.

  <p>The problems described before for boundary faces intersected by a
grid edge we have also for boundary lines intersected by grid faces.
This case is more complicate to handle, because it will be not so easy
to detect or avoid it in the previous steps, it is also not so easy to
detect it in the Delaunay grid, and it is often more complicate to
rebuild.