Boundary conditions are specified on <a href="geom.html#internal-boundaries">internal boundaries </a> (both 
multiple-region and
repeated-boundary
 internal boundaries) in the same way as on ordinary boundaries.
An internal boundary with a Dirichlet condition will force <em>u</em> to take
on the prespecified value on that boundary, and hence <em>c&#183;du/dn</em>
 will
probably be discontinuous.

<p>
A Neumann condition of 0 on an internal boundary (the default) means
that the internal boundary does not act as a boundary at all, i.e.,
the nodes on that boundary are treated like interior nodes.

<p>
A nonzero Neumann condition on an internal boundary is treated like
a prespecified jump in <em>c&#183;du/dn</em>, whereas <em>u</em>
 remains continuous across
the boundary. 

<p>
<a name = "double"></a>
Sometimes it is desirable to have separate boundary conditions
on the two sides of an internal boundary.  For this purpose,
a procedure for doubling internal boundaries is available:
<code>gmdouble</code>.  It has two calling sequences:

<blockquote>
Matlab:
<strong>
<em>newbrep</em> = gmdouble(<em>brep, facelist</em>);<br>
</strong>
or:
<strong>
[<em>newbrep, newmesh</em>] = gmdouble(<em>brep, facelist, mesh</em>);<br>
</strong>
Tcl/Tk:
<strong>
gmset <em>newbrep</em> [gmdouble <em>brep facelist</em>]<br>
</strong>
or:
<strong>
gmset {<em>newbrep newmesh</em>} [gmdouble <em>brep 
facelist mesh</em>];<br>
</strong>
</blockquote>


This routine
takes a brep whose regions have internal boundaries.
Some subset of those internal boundaries 
is specified by <em>facelist</em> by name;
this is a list of strings specifying the names of 
faces of dimension <em>d</em>&#8722;1 
that are internal boundaries to be doubled.
In matlab, this argument is a cell-array of strings.  In Tcl/Tk,
it is a list of strings.

<p>
These must be repeated-boundary internal boundaries.
The routine replaces all the specified internal boundaries
with two copies.  It may also duplicate lower dimensional
faces that interconnect the specified internal boundaries.  The
resulting brep is the return argument <em>newbrep</em>.

<p>
The <em>facelist</em> array can optionally be an asterisk (specified
as {'*'} in matlab, since the argument must still be a cell
array of strings).  In this case, the routine will duplicate
all repeated-boundary internal boundaries.

<p>
The second calling format takes a mesh as the third argument.
This argument is a mesh generated for the
original brep by <code>gmmeshgen</code>.  The <code>gmdouble</code>
routine returns
another mesh in which mesh nodes are duplicated
if they sit on faces of the brep
that were duplicated.  The simplices in the mesh are
adjusted accordingly.  The resulting mesh
is returned in argument <em>newmesh</em>.

<p>
You can now attach different boundary conditions to the two copiesof the doubled faces.  To help
you figure out which of the two doubled faces points in which direction,
you can use the orientation information of the faces.

<p>
Note that you should not run the mesh generator on
the brep produced by <code>gmdouble</code>
because it is an invalid brep (there are two coincident faces).
But it is OK&#8212;in fact, recommended&#8212;that you run <code>gmchecktri</code>
 on the new brep-mesh pair.
Thus, for a domain with doubled internal boundaries,first create them as single-layer internal boundaries, then
generate a mesh, then execute <code>gmdouble</code>, then 
execute <code>gmchecktri</code>, and then finally
execute <code>gmfem</code>.


<p>
The nodes and simplices of the new mesh are in correspondence
with the nodes and simplices of the previous mesh according to
the following rules.  The order of simplices associated with
each face is preserved.  The nodes that are duplicated occur
together.  Duplicates of (global) node number <em>i</em> in thenew mesh are numbered <em>i</em>, <em>i</em>+(<em>n</em>+1),
<em>i</em>+2(<em>n</em>+1), etc., where <em>n</em> is the
maximum global node number in the original mesh.  This 
system prevents any numbering conflicts.

<p>
The <code>gmdouble</code> routine uses a tolerance for casting
rays to determine which side of the internal boundary contains
mesh points.  There is no way in the current calling sequence
to specify the tolerance, so <code>gmdouble</code> always
uses the <a href="meshgen.html#default_tol">default tolerance</a>.

<h2> Numerical methods used </h2>

The BVP is discretized using piecewise linear
basis functions.  The Dirichlet data is linearly interpolated from
nodal values, and the Neumann data is integrated using a low-order
quadrature rule that evaluates the Neumann data at mesh node points
only. The source data is integrated using a low-order quadrature
rule that evaluates the source function at interior and boundary nodes.
The conductivity is integrated on interior elements using
the midpoint rule, and hence the conductivity function is evaluated
only in the interior.

<p>
For a BVP in which <em>c</em> is constant, <em>f</em>=0, the
solution <em>u</em> is linear, and the domain is
polyhedral, the approximations made in the
last paragraph are exact, so the package should return
the exact value of <em>u</em> up to rounding error for any combination
Dirichlet and Neumann boundary conditions.

<p>
The stiffness matrix and Dirichlet boundary data are assembled by
forming triple products <em>K=A<sup>T</sup>DA</em> and 
<em>f<sub>0</sub>=A<sup>T</sup>Db</em>, where <em>A</em> is matrix that
depends only on the mesh, and <em>D</em>, 
a diagonal matrix, depends on the
conductivity field.  The vector <em>b</em> is depends on the boundary data,
source term, and the mesh.

<p>
The matrix <em>K</em> is symmetric, positive definite, and sparse.  The
number of rows and columns of <em>K</em> is equal to the number of nodes
in the mesh that are not Dirichlet boundary nodes.

<p>
This method for assembling the matrix is nonstandard but is
mathematically equivalent to the standard method.  
(The standard method involves forming ``element stiffness matrices''
and then adding them to get <em>K</em>.) The method is
described in detail in section 3 of 
<blockquote> 
S. Vavasis, ``Stable
finite elements for problems with wild coefficients,''  SIAM
J. Num. Anal.  33 (1996) 890.
</blockquote> 
It should be noted that the main result of
that paper, which is an algorithm for accurately solving the finite element
linear system <em>Ku=f</em> for BVP's in which <em>c</em> varies by orders of magnitude
over the domain,  has <em>not</em> been implemented in the current
package.

<p>
The advantage of assembling the matrix <em>K</em> as 
<em>A<sup>T</sup>DA</em> is that
all the operations vectorize well in Matlab.  There are no explicit
loops over elements or nodes anywhere in the Matlab
finite element m-files
so it performs well.

<p>
The linear system <em>Ku=f</em> is solved via the Matlab backslash command
applied to the sparse matrix <em>K</em>.  By default this is sparse symmetric
minimum degree ordering.  

<p>
The ordering used by <code>gmfem</code> by default is the ordering that emerges
from the mesh generator together with the setting of <code>spparms</code> 
(which
defaults to minimum degree).


<hr>
<p>
This documentation is written by 
<a href="http://www.cs.cornell.edu/home/vavasis/vavasis.html">Stephen A. 
Vavasis</a> and is
copyright &#169;1999 by 
<a href="http://www.info.cornell.edu/CUHomePage.html">Cornell 
University</a>.
Permission to reproduce this documentation is granted provided this
notice remains attached.  There is no warranty of any kind on
this software or its documentation.  See the accompanying file
<a href="copyright.html">'copyright'</a>
for a full statement of the copyright.
<p>
<address>
Stephen A. Vavasis, Computer Science Department, Cornell University,
Ithaca, NY 14853, vavasis@cs.cornell.edu
</address>
</body>
</html>