thgdalg.html

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

<h1><a NAME="gdGReal">Algorithms for Contravariant Geometry Descriptions </a></h1>


In this part we consider different algorithms which allow to define
cogeometries.

 
<h2><a NAME="subdiv">Simplex Subdivision</a></h2>


We have already found that the cogeometry may be defined by their
results for arbitrarily small simplices. So, it makes sense to consider
also algorithms which do not work correctly for big simplices.
Assume we have such an algorithm and some epsilon so that the
algorithm may be considered as correct for smaller simplices.

  <p>The subdivision algorithm allows to obtain the correct result also
for greater simplices. It works so:

 
 <ul>
 <li>If the simplex is smaller than epsilon, the given algorithm will be called.
 
 <li>Else, the k-simplex will be subdivided into 2^k sub-simplices of the
half size.
 
 <li>It must be detected which simplex contains the initial flag. In
this step, degenerate cases have to be handled if the initial flag is
in a neighbourhood of a boundary between sub-simplices.
 
 <li>Then a cycle over the sub-simplices has to be considered. For the
given sub-simplex, the same algorithm will be called recursively. If
the result is a k-flag inside the sub-simplex or a (k-1)-flag on the
outer boundary, this result will be returned. Else, the result is a
(k-1)-flag on a boundary between two sub-simplices. In this case, in
the next step we consider the related neighbour sub-simplex and use the
result flag of the previous step as input.
 
 </ul>


Is it possible to get an infinite loop in this algorithm? In the
general case, there will be only a finite number of intersections of
the related boundaries with the inner boundaries between the
sub-simplices. So, an infinite loop may be only a cycle between a
finite set of flags. If the initial algorithm is really symmetric,
this is not possible. Thus, usually there will be no infinite loops,
but degenerated situations and errors of the initial algorithm may
lead to such infinite loops.

  <p>In a typical regular situation, the number of intersections with
inner boundaries will be small. Thus, if the number of steps will be
great, an infinite loop may be assumed. Thus, it seems useful to break
the loop after a maximal number of steps which may be not very
large. The last (k-1)-flag after a break lies inside the simplex. The
value may be returned as a k-flag on some artificial "error boundary".
This allows to continue the computation.

  <p>Another idea for the error handling is to restart the computation
with a temporary smaller value of epsilon. This seems useful, because an
incorrect value for epsilon seems to be the most probable error.

 
<h2><a NAME="default">The Default Function</a></h2>


Assume we have defined the first (k-1) functions f(i). Then we can
define a default function f(k). This function subdivides the
(k-1)-boundary into parts with identical higher-dimensional
neighbourhood. The resulting k-boundary consists of a unique default
k-segment.

 
 <ul>
 <li>We consider the border of the simplex and start with the given
(k-1)-flag on the first side. We consider the continuation of this
flag on the side of the simplex.
 
 <li>While we obtain a (k-2)-flag on the border of a side, we continue the
search on the neighbour side.
 
 <li>Else, we have found another (k-1)-flag on some side of the
simplex. If this flag is part of the same (k-1)-segment, it will be returned as the continuation of the initial flag through the simplex.
 
 <li>Else, we assume that there is some point of a k-segment inside
the simplex which subdivides the (k-1)-boundary line inside the
simplex into two different (k-1)-segments which ends we have found on
the simplex border.
	
 <ul>
 <li>If the size of the simplex is less than the required
accuracy epsilon we return a k-flag which consists of the middle point
of the simplex as the approximation of a point of the default
k-segment and the given (k-1)-flag on the side.
        
 <li>Else, we make a further recursive subdivision of the
simplex in analogy to the previous algorithm.
	
 </ul>

 
 </ul>


Possible infinite loops for the search over the simplex sides may be handled
in analogy to the previous algorithm.

  <p>In the case of multiple intersections of the (k-1)-segment with
the border of the simplex, it is possible that the algorithm finds not
the first intersection of this segment with the border. Such a result
is incorrect. But this error usually occurs only for big
simplices. That's why further subdivision may be used to avoid such
errors.

  <p>This default function f(k) returns only one unique "default
segment", it does not allow to transfer nontrivial attributes of the
k-boundary. The default function f(k+1) will not lead to any
subdivision of this default k-segment, that means it does not create
any (k+1)-flag. So, there will be no possible input values for the
function f(k+2), the geometry has the codimension k. Thus, the
implementation of the default functions f(k) and f(k+1) allows to
use the first functions f(i) for i &lt; k of a cogeometry as a
correct, complete cogeometry of codimension k.

  <p>This leads to a powerful "fast prototyping" strategy. In the first
step, only the function f(0) has to be implemented. Later, the other
functions f(k) may be implemented step by step. Every step leads to
a better, more accurate and more powerful description of the resulting
geometry, but even the first step leads to a complete geometry
description which may be used as a prototype of the correct
cogeometry.

  <p>The algorithm may be slightly modified for the case of f(2).
Even if there is no subdivision of the boundary face into
parts, boundary faces with different pairs of neighbour regions may be
considered as different. So, a geometry with nontrivial boundary lines
(that means with codimension two) may be created even if only f(0) is
defined.

  <p>The other flag points in the algorithm are in the simplex. That
means, they usually will be not orthogonal to the boundary.

  This algorithm was one of the reasons to use flags instead of
intersection points as in the first variant and the implementation in
IBGD. The (k-1)-part of the initial flag is necessary in this
algorithm.
 

 
<h2><a NAME="induced">The Induced Cogeometry</a></h2>


As we have already mentioned, if we have a smooth mapping f: X \to
Y and a cogeometry G(Y) on Y, we can define on X an induced cogeometry
G(X). We have remarked, that it is a failure of the standard
geometry description that there is no general algorithm which allows
to create this induced geometry. Let's consider now how the cogeometry
allows to describe the induced geometry.

Let's consider at first the case of an affine mapping f between affine
spaces. We have to define the function f(k) on X for given functions
f(i) on Y. For a given input flag and the related side or simplex,
their image also defines a correct, non-degenerated pair flag - side
or simplex, because every flag obtained by the following algorithm
will be created so that this property is fulfilled. We compute the
resulting flag on Y. The pre-image of the resulting point intersects
the related side or the initial simplex in a single point, which will
be used as the position of the resulting flag. The other flag points
may be defined in the same way (nonorthogonal variant) or so that the
related directions d(i) are orthogonal to the pre-image of the
position of the flag in Y (orthogonal variant).

  <p>To compute the related positions it is in the nonorthogonal case
not necessary to compute something about the pre-image. This may be
dangerous if a rounding error moves the point out of the image. We can
simply use the barycentric coordinates of the resulting points in the
image simplex in Y to compute their position in the original simplex
in X.

  <p>Thus, we have an easy, straightforward algorithm for the affine
case. For the smooth case, we simply can use subdivision into smaller
simplices until the simplices are small enough to allow the
approximation by the affine algorithm.

  <p>We can modify the algorithm to make it faster: Instead of
subdividing the simplex into equal parts, we can use the barycentric
coordinates of the result flag in Y to subdivide the simplex.

  <p>A variant of the algorithm may be used for piecewise affine
mappings. We explicitly subdivide the simplex into pieces so that the
mapping is affine on every piece, and use the affine variant on each
piece. This seems to be useful for mappings with a small number of
different parts, for example for the continuation of a geometry
defined in a cube to the outside.

 
<h2><a NAME="intersect">The Intersection of Geometries </a></h2>


The intersection of geometries is another example of a natural
operation which is hard to implement for the standard geometry
description. let's look how the intersection may be implemented in the
contravariant geometry description.

  <p>In general, a k-segment S_k of the intersection will be
described by an i-segment S_i of the first and an (k-i)-segment
S_k-i of the other cogeometry. The codimension of the resulting
geometry is the sum of the codimensions of the two cogeometries.

  <p>The general scheme will be analogical to the case of an induced
geometry. We can use subdivision until the simplices are small enough
to be approximated by the affine situation. To consider the affine
situation for some fixed pair (k, i) is straightforward but complicate
to implement.

  <p>For the first functions the implementation is much more trivial:

 
 <ul>
 <li>f(0) is very simple. We have to get the two regions containing
the point in the two cogeometries, and the resulting pair defines the
region in the intersection.
 
 <li>f(1) is also simple to implement. We call f(1) for above
geometries and use the first intersection.