thgdcov.html

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

 <li>The position of the flag is the position of k.
 
 </ul>


There are some advantages of this method:

 
 <ul>
 <li>We need only one data type for the description of all flags.
 
 <li>If we need only the intersection point instead of the complete
flag, we can simply ignore the information which is not necessary.
 
 <li>Sometimes it is useful to have a point inside a segment in the
immediate neighbourhood of the boundary. But even a small shift from
the basic point into the related direction can lead to a point of
another segment. So, to obtain such a neighbour point becomes easier.
 
 <li>The point data type may be used to transfer function values
through the geometry description.  If the function we have to consider
has different boundary limits for different segments, a point inside a
given segment but in the neighbourhood of the boundary point may be
used to transfer the boundary limit.
 
 </ul>


 
<h3><a NAME="rounderror">Rounding Error Handling</a></h3>


In an implementation of the interface we have to manage two problems:

 
 <ul>
 <li>degenerate cases and
 
 <li>rounding errors.
 
 </ul>


  <p>The main problem is that the result of a call may be used for
later input. So, the result must have properties which allow the
continuation of the computation if it will be used later as input.
The typical situation is that we have some expression f which in the
continuous situation leads to one answer if f &gt; 0 and to another if
f &lt; 0, in the case f = 0 the answer is not defined. In the case of
finite precision, not only the case f = 0 is a problem, but also very
small values of f, because different rounding errors for f may lead to
a different classification of the situation in different parts of the
algorithm. This usually leads to a fatal error in the program.

  <p>The strategy to avoid such situations is to make a small
modification of the result if the exact result is in such a dangerous
neighbourhood of a degeneration. The modification must be small enough
compared with the required accuracy of boundary computation, but it
must be large enough to avoid an incorrect classification if it will be
used later as input.

  <p>Thus, if the required accuracy is large enough compared with the
possible rounding errors, this technique allows to avoid fatal
errors. It also does not require a special handling for the
degenerated situations where the result is not defined in the exact,
continuous case.

  <p>Remark that the case of a degenerated simplex is not dangerous, if
the side containing the input flag itself is not degenerated. There
will be simply a smaller set of possible output --- there may be no
k-flag inside the simplex and no (k-1)-flag on degenerated sides.

 
<h3><a NAME="cpptwofunc">Subdivision into Two Different Functions</a></h3>


For a theoretical consideration it looks very nice if we have only one
function for every dimension. But in the real implementation it
becomes easier two distinguish two functions:

 
 <ul>
 <li>The first variant of f(k) for the case of a (k-1)-flag as input.
 
 <li>The other variant of f(k) for the case of a k-flag as input.
 
 </ul>


The idea of the simplification is that for the first variant we
implement only one special case --- the flag on the first side of the
simplex. This makes the implementation simpler, but for the calling
function it is not difficult to use a point order for the simplex
corner so that the flag lies on the first side. For the other variant,
we can use a default implementation:

  Subdivide the simplex into smaller simplices so that the flag
lies on the border between the sub-simplices. Use the first variant of
f(k) to find the continuation. While the continuation was found on
the inner border between the two sub-simplices, we have to continue the
search in the other sub-simplex.
 

 
<h3><a NAME="nonorthflag">Nonorthogonal Flag Directions</a></h3>


The orthogonality condition for the flag directions require a special
consideration.

 
 <ul>
 <li> The orthogonality may not be exactly fulfilled in finite
precision arithmetics caused by rounding errors. Thus, in reality we
will not have exactly orthogonal flag directions.
 
 <li> There are algorithms which do not include the computation of
the related orthogonal directions. Usually they allow to create only some
set of directions with non-degenerated projection on the simplex plane.
 
 <li> For a non-smooth boundary, it will not be possible to define the
tangential and orthogonal directions required for the definition of a
flag.
 
 </ul>


These problems may be solved using the convention that the directions
must not be orthogonal, but only their projections have to be not
degenerated. But in this case we obtain a new problem --- the
projection of the directions on the simplex plane may lead to an
incorrect result for the orientation of the intersection. This problem
may be solved by the following convention:

 
 <ul>
 <li>The flag directions of the result must define a set of
non-degenerated projection on the simplex plane so that it's
orientation coincides with the projections of the orthogonal flag
directions.
 
 <li>For every flag, a plane containing this flag is defined by the
following rules:
	
 <ul>
 <li>If the flag is immediate output of some function call,
the plane is defined by the simplex or it's side containing the flag.
	
 <li>If the flag l is a part of another flag k with k
&gt; l, the plane is the part of the plane of the initial flag which
is orthogonal to the flag directions d(k),...,d(l+1) which will be
omitted.
        
 </ul>

 If the flag will be used as input, the plane of the simplex
containing this flag must coincide with this plane containing the
flag.
 
 </ul>


Very often the flag will be used only in combination with the same
simplex. The previous rules allow to use also other simplices if these
simplices are in the same plane. Our general strategy to handle
rounding errors guarantees that the orientation of the projection will
be the same also for a slightly different plane caused by rounding
errors.

  <p>We see that the orthogonal variant is easier to use because it
allows to use any plane which contains the given flag and leads to a
non-degenerated projection of the flag direction. It is not necessary
to hold information about the plane containing the given flag.

 
<h3><a NAME="cppinterf">C++ Interface and OOP</a></h3>


The contravariant geometry description is very natural from
object-oriented point of view:

 
 <ul>
 <li>The cogeometry is a class.
 
 <li>The intersection functions f(k) are the methods of this class.
 
 </ul>


This leads to a natural C++ implementation. The cogeometry will be
defined as an abstract class, the functions f(k) will be virtual
methods.  A concrete geometry description will be a derived class
which contains all necessary data of the geometry. In principle, for
all functions f(k) there is some default implementation which may be
used if no better implementation is available. This default
implementation describes a boundary which is not further subdivided
into parts.

 
<h3><a NAME="cppsimpl">Special Conventions for Codimension 0 and 1</a></h3>


The case f(0) may be in principle considered as the degenerated case
of the general intersection function f(k) for k = 0, but this is not
very useful because there is no (-1)-flag and there are no sides for
the 0-simplex. So, let's define f(0) separately:

 
 <ul>
 <li>For a given point (a 0-simplex) f(0) returns a 0-flag on this point.
 
 </ul>


That means, it returns a pair which consists of the initial point and
the region containing this point. The nontrivial part of this operation is
to find the region containing the given point.

  <p>The function f(1) in principle may be used as defined, but we
use here also some modification. Every smooth boundary face lies
between two regions. Thus, if for every higher codimension the number
of flags containing a given boundary point is not bounded, it is
well-defined for codimension 1. That's why it seems natural to return
above flags in a single call of f(1) instead of calling f(1)
always two times from above sides of the edge. We simply have to add one
new output point --- the point 0 of the "other side" of the flag.

<h2><a NAME="gdattrib">Attribute Handling</a></h2>


Necessary part of the description of a physical situation is the
attribute description. An <b>attribute</b> may be an arbitrary
application-dependent information. An attribute has a result type
(scalar, vector, integer) and can depend on geometrical data (points,
segments). The physical sense of the attribute is defined by the
application and hidden from the geometry description. That's why it is
a good idea to separate the geometry description and the attribute
description whenever possible. For many types of attributes this is
possible.

  <p>Consider, for example, attributes which depend on the segment but
are constant on this segment like material data of a region or
boundary conditions on a boundary face. In this case, the geometry
description returns a pointer (or a small positive integer) as the
identifier of the segment, and the related attributes may be described
by components of the related object (or by arrays over the possible
integer values). In this case, we have a clear modular subdivision
between geometry description and application data.

  <p>Another interesting class of attributes are space-dependent
attributes. In principle, they may be described independently of the
geometry description. But often this leads to problems. The reason is
that they have to be evaluated once for every point which will be
considered in the application. If this may be done later, this does
not cause problems. But often the values may be used in the process
itself. For example, the refinement criterion used for grid generation
may depend on such an attribute. That's why the attribute must be
evaluated in the process of grid generation. Or the attribute value of
an initial geometry may be used to define another geometry (for
example the induced geometry of a mapping defined by this attribute).
In this case, the attribute must be evaluated in the process of
geometry description itself. Thus, a separate description may lead to
problems.

  <p>This situation is a special case of another class of attributes
where it is not so easy to separate the geometry description and the
attribute description. These are attributes which depend not only on
the segment, but also on the point of the segment. This class of
application also often occurs in applications. For example, many
functions are defined only on some regions. Boundary concentrations
and other functions defined only on a boundary also often have to be
considered. Another example are functions which are continuous only
inside a given region, but with different boundary limits for
different regions (for example concentrations with segregation
boundary conditions).

  <p>Because of this deep interaction with the geometry description it
is necessary to have a general scheme to manage such attributes.
Let's consider now this interaction. At first, remark that we have
also a natural functional behaviour for the attributes. Indeed, for
the geometry induced by a mapping f:X --> Y on X and for given
attributes of the original geometry on Y related induced attributes
may be defined in a natural way: The attribute value of a point in X
is simply the attribute value of it's image in Y.

  <p>For the standard geometry description, such attribute data may be
described by functions over the cells of the cell complex. The
resulting attribute description is not very modular --- a small
modification of the mapping (another parametrization) requires a
modification of the related attribute description on the cell. We also
obtain the same incorrect functional behaviour as for the cell complex
itself. Indeed, a mapping f: X --> Y and an attribute description
on a cell complex in X defines a natural attribute description on the
image of the cell complex in Y. So, not only the construction of the
cell complex which describes the induced geometry, but also the
transfer of the attribute description becomes a complicate problem.

  <p>In the case of the contravariant geometry description, we can manage
attributes using the following simple technique:

 We use a data type for the point which contains also the
attribute values of the point. The function f(k) has to compute also
these attribute values for the output points.
 

  <p>In principle, this technique may be considered as a special case
of an induced geometry: We have to consider the embedding of a
lower-dimensional space into a higher-dimensional space defined by the
graph of the attribute functions.  This interpretation shows that it
is possible to combine this technique with other operations using the
composition of mappings. For example, the attribute values of some
geometry may be used to define a mapping which may be used to define
an induced geometry.

  <p>Our scheme leads to a natural implementation for the
interpolation of function values for a given grid. The functionsf(k) have to find the intersections of a simplex with the related
boundary.  Using our technique, this function also has to evaluate the
function values.  But this is a very natural place to interpolate the
function values, because the majority of data we need for the
interpolation of the function values in the grid (especially the
element which contains the point) we need also if we have to find only
the intersection point in the grid.