dolfin team mailing list archive

["Garth N. Wells" <gnw20@xxxxxxxxx>]

On 28/03/11 12:11, Anders Logg wrote:
> On Mon, Mar 28, 2011 at 10:37:10AM +0100, Garth N. Wells wrote:
>>
>>
>> On 25/03/11 12:58, Anders Logg wrote:
>>> On Fri, Mar 25, 2011 at 10:28:20AM +0000, Garth N. Wells wrote:
>>>>
>>>>
>>>> On 24/03/11 19:10, walker@xxxxxxxxxxxx wrote:
>>>>> On Thu, March 24, 2011 1:32 pm, Anders Logg wrote:
>>>>>> On Thu, Mar 24, 2011 at 05:19:03PM +0000, Garth N. Wells wrote:
>>>>>>>
>>>>>>>
>>>>>>> On 24/03/11 14:28, walker@xxxxxxxxxxxx wrote:
>>>>>>>> On Thu, March 24, 2011 8:53 am, Peter Brune wrote:
>>>>>>>>> Ah, this argument, again.  I continue to disagree.  While you could
>>>>>>>>> certainly home-roll a FEM code that does this under assumptions like
>>>>>>>>> constant quadrature order and small meshes,  this is going to fall
>>>>>>> down
>>>>>>>>> hard
>>>>>>>>> as soon as you have varying quadrature degree or type, which we
>>>>>>> always do.
>>>>>>>>> This would also keep around way more state than the current assembly
>>>>>>>>> paradigm allows, not to mention the per-quadrature-point storage
>>>>>>>>> requirements of a full-blown tensor.
>>>>>>>>
>>>>>>>> I don't see why constant quadrature order is so terrible.  Is it
>>>>>>> really
>>>>>>>> necessary to have varying quadrature degree?  What do you mean be more
>>>>>>>> state info?  This is just another locally defined function.  I am a
>>>>>>> little
>>>>>>>> unclear about what you are saying.
>>>>>>>>
>>>>>>>
>>>>>>> I don't see what the 'argument' is either. What is it exactly?
>>>>>>>
>>>>>>>>> The other problem with this approach is that you severely restrict
>>>>>>> the
>>>>>>>>> spaces your coordinates can live in by having some oddly-defined
>>>>>>> notion of
>>>>>>>>> "higher-order vertices.  If you have some coefficient that is the
>>>>>>>>> coordinates instead it's automatic.  I think you definitely should
>>>>>>>>> recreate
>>>>>>>>> the jacobian in tabulate_tensor; it's what we do now; it's what we
>>>>>>> should
>>>>>>>>> do
>>>>>>>>> if we don't want to store per-quadrature-point tensors on the DOLFIN
>>>>>>> side
>>>>>>>>> when the DOLFIN side doesn't even know about quadrature points.
>>>>>>>>
>>>>>>>> I am not proposing having higher order vertices.  I am saying treat
>>>>>>> the
>>>>>>>> geometric map like **it has its own finite element space**!  If you
>>>>>>> want
>>>>>>>> to have a special map for each bilinear form (so you can have varying
>>>>>>> quad
>>>>>>>> degree) then fine.  But keep in mind that once you have a non-linear
>>>>>>> map
>>>>>>>> for the geometry, you will (in general) not be able to compute the
>>>>>>> tensors
>>>>>>>> exactly.  So you will need to fix the quad degree somehow.
>>>>>>>>
>>>>>>>
>>>>>>> I agree with this. Forgetting forms and PDEs, we should be able able to
>>>>>>> have a pure geometric representation. Say, for example, that I have a
>>>>>>> mesh that uses some mapping and I want to know (possibly approximately)
>>>>>>> its volume. The mesh and cells need to know all the geometry. Also, say
>>>>>>> I want to know in which cell a point (in the real coordiantes) lies.
>>>>>>> This is not related to forms.
>>>>>>
>>>>>> We've had this argument many times before... The thing that I've
>>>>>> always found attractive with Peter's approach is that one can
>>>>>> represent the geometry (the embedding of the mesh) as a field on top of
>>>>>> a standard mesh (with straight edges). So it doesn't require adding
>>>>>> tons of new data to the mesh class.
>>>>>>
>>>>
>>>> How would this work for very complicated meshes, e.g. a car, engine,
>>>> etc? In these cases one doesn't have a straightforward functional
>>>> description of the geometry, but just points that lie of the boundary or
>>>> a complex collection of spline.
>>>
>>> We would then need a function space that can represent those kinds of
>>> fields (NURBS). If we add such functionality (to DOLFIN and FFC), we
>>> could also use the NURBS as basis functions ("isogeometric analysis")
>>> which seems to be popular these days.
>>>
>>
>> This is missing the point. Complex geometry is given, and not
>> constructed for convenience by the person performing the analysis. It
>> can come in many formats. It is necessary to support geometries that are
>> provided.
> 
> No one is arguing against that. The discussion is how we should are lotsstore
> the geometry: either in a special MeshGeometry class or as a Function.
>

This is an over simplified perspective. We are not discussing whether or
not to use Function, which since we all agree of a functional
representation is an implementation detail. The question is what is
required to represent a broad range of geometries.  Splines (the types
of interest) are not defined on simplices, so it's not a simple
questions of straightforward extension of FIAT. To evaluate splines,
geometric data is required.

(I'm wishing that NURBS were never mentioned ;). There is more out there
than just NURBS - we need to deal with whatever is available. Part of my
assertion is that FFC + FIAT cannot possibly support everything.)

> My point is that if we store it as a Function (which among other
> things requires adding the corresponding basis functions in FIAT) then
> we could not only represent the geometry using NURBS, but we could
> also use them to compute. Another advantage would be that we wouldn't
> need to extend the MeshGeometry class or create a new MeshGeometry
> class to handle the new geometry.
>

How could we construct NURBS on a domain without some geometric data for
constructing the NURBS? It's like an FE simulation on a domain, but
without a mesh to define the domain.

Garth


>> Obviously, all representations must have some form of
>> functional representation, otherwise it's not possible to compute the
>> Jacobian at arbitrary points. Having a functional representation is
>> something that we all agree on, but the interface and code generation
>> approach must recognise the complexity and diversity of geometric
>> models, and that we don't get to choose the geometric model.
> 
> No one is arguing against that either.
> 
> --
> Anders
> 
> 
>> Garth
>>
>>>>>> I've been attracted to this approach ever since Matt first presented
>>>>>> the design for the Sieve where one of the main points is to representl
>>>>>> the coordinates of the mesh as a field, or in Sieve-speak a section of
>>>>>> a vector bundle... ;-)
>>>>>>
>>>>>> But this argument may be worth less these days as we can no longer
>>>>>> store fields (Functions) to file.
>>>>>>
>>>>>
>>>>> Is this really different from having a finite element space for the
>>>>> mapping?  You can define a finite element function in that space to store
>>>>> the higher order mesh data.  Or maybe have a separate higher order
>>>>> geometry class that can do that OR some other representation like NURBS.
>>>>>
>>>>
>>>> This is what I would like - the flexibility to use a functional
>>>> representation (including functions other than FE functions, such as
>>>> NURBS) or an interpolated approach with higher-order vertices.
>>>>
>>>> Garth
>>>>
>>>>> - Shawn
>>>>>
>>>>>>> Garth
>>>>>>>
>>>>>>>
>>>>>>>>> Arguments about the interface can come later.  I'm sanity checking
>>>>>>> the
>>>>>>>>> code
>>>>>>>>> I updated to the latest FFC with a couple test cases and should have
>>>>>>> it in
>>>>>>>>> a
>>>>>>>>> repo somewhere public shortly.
>>>>>>>>>
>>>>>>>>> - Peter
>>>>>>>>
>>>>>>>> I think I will let other people figure out how to put this into
>>>>>>> FEniCS.
>>>>>>>> Obviously, my experience is limited.  However, my code is available on
>>>>>>> my
>>>>>>>> web page: https://www.math.lsu.edu/~walker/FELICITY.html
>>>>>>>>
>>>>>>>> I am NOT advertising my code.  I am just offering it as an example of
>>>>>>> what
>>>>>>>> to do or what NOT to do.
>>>>>>>>
>>>>>>>> If you are curious about how the generated code looks, please have a
>>>>>>> look.
>>>>>>>>  All you need to have is MATLAB and a C++ compiler.  So the
>>>>>>> instructions
>>>>>>>> are:
>>>>>>>>
>>>>>>>> 1. Open MATLAB.
>>>>>>>> 2. type "mex -setup" to pick the C++ compiler (this is easy).
>>>>>>>> 3. unzip my code to a directory.
>>>>>>>> 4. within MATLAB, change to that directory.
>>>>>>>> 5. run this script: "FELICITY_paths.m".  This will add all of the
>>>>>>>> appropriate paths and unit tests.
>>>>>>>> 6. then run this script: "test_FELICITY.m".  This runs all the unit
>>>>>>> tests
>>>>>>>> I have so far.
>>>>>>>> 7. It you want to see the code that is generated, then just go to the
>>>>>>>> Unit_Test subdir under the Matrix_Assembly dir.
>>>>>>>>
>>>>>>>> Steps 1-6 should not take more 5 mins (maximum).  I have tested it on
>>>>>>>> windows and LINUX.  If it doesn't run, please let me know.
>>>>>>>>
>>>>>>>> - Shawn
>>>>>>>>
>>>>>>>>> On Thu, Mar 24, 2011 at 7:09 AM, Garth N. Wells <gnw20@xxxxxxxxx>
>>>>>>> wrote:
>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> On 24/03/11 04:38, walker@xxxxxxxxxxxx wrote:
>>>>>>>>>>> Hello Anders and Garth.
>>>>>>>>>>>
>>>>>>>>>>> You may not remember me, but I had some interest way back in seeing
>>>>>>>>>> this
>>>>>>>>>>> higher order geometry stuff put in.
>>>>>>>>>>
>>>>>>>>>> I remember!
>>>>>>>>>>
>>>>>>>>>>> However, it seemed like P. Brune had
>>>>>>>>>>> a way so I backed off.   In the meantime, I decided to stick with
>>>>>>> my
>>>>>>>>>> own
>>>>>>>>>>> FEM toolbox setup and continued developing it.  I have higher order
>>>>>>>>>>> elements implemented (in principle for any finite element that is
>>>>>>>>>> defined)
>>>>>>>>>>> and I can do this for 1-D curves in 3-D and 2-D surface meshes in
>>>>>>> 3-D
>>>>>>>>>>> (which is another thing I would like to see in FEniCS), in addition
>>>>>>> to
>>>>>>>>>> 2-D
>>>>>>>>>>> and 3-D curved meshes.  I can even compute things like the total
>>>>>>>>>> curvature
>>>>>>>>>>> vector (i.e. 2nd fundamental form).  I even implemented some code
>>>>>>>>>>> generation in my setup because I think you have a good idea there.
>>>>>>>>>> But
>>>>>>>>>> my
>>>>>>>>>>> goal is not to duplicate FEniCS (there is no way I could do that).
>>>>>>> I
>>>>>>>>>> just
>>>>>>>>>>> needed these basic components for my own work.
>>>>>>>>>>>
>>>>>>>>>>> Anyway, I would like to describe how I did it, b/c the lessons I
>>>>>>>>>> learned
>>>>>>>>>>> may be useful to you.
>>>>>>>>>>>
>>>>>>>>>>> 1. Garth is correct.  All local geometric transformations should be
>>>>>>>>>>> computed separately from the tabulate tensor routine.  It has been
>>>>>>>>>> awhile
>>>>>>>>>>> since I perused the FEniCS code, but I assume that is still the
>>>>>>> same.
>>>>>>>>>> I
>>>>>>>>>>> repeat, you should NOT compute the jacobian (or derivatives of the
>>>>>>>>>>> jacobian, or the normal vector, etc...) in the tabulate tensor.
>>>>>>> This
>>>>>>>>>> has
>>>>>>>>>>> the following advantages: (a) the code is more modular; you
>>>>>>> separate
>>>>>>>>>> the
>>>>>>>>>>> geometry transformation from the element calculation; (b) you can
>>>>>>>>>> *reuse*
>>>>>>>>>>> the jacobian calculation for each bilinear form that has to get
>>>>>>>>>> assembled,
>>>>>>>>>>> instead of recomputing it every time you enter tabulate_tensor.
>>>>>>>>>>>
>>>>>>>>>>> I personally found this to work great for me.  And I would like to
>>>>>>>>>>> emphasize that I do not see a 10-fold slow-down in assembly time
>>>>>>> (for
>>>>>>>>>> a
>>>>>>>>>>> fixed quad rule); I don't remember the timings, but it seemed to be
>>>>>>>>>>> negligible b/c this calculation is done in cache.  Again, the local
>>>>>>>>>>> mapping transformation for each element should have *nothing to do*
>>>>>>>>>> with
>>>>>>>>>>> the bilinear forms, etc.  It is a separate calculation (of course,
>>>>>>> you
>>>>>>>>>>> need to take into account Piola transform when appropriate, but
>>>>>>> this
>>>>>>>>>> is
>>>>>>>>>>> minor).
>>>>>>>>>>>
>>>>>>>>>>> 2. When the user wants to specify the type of "higher order"
>>>>>>> geometry,
>>>>>>>>>> I
>>>>>>>>>>> would suggest the following.  You should have the user define a
>>>>>>>>>>> "Geometric_Element" just like when they define regular finite
>>>>>>>>>> elements,
>>>>>>>>>>> i.e.
>>>>>>>>>>>
>>>>>>>>>>> Geometric_Element = FiniteElement("triangle",1,1)
>>>>>>>>>>>
>>>>>>>>>>> Or something like that (I don't remember your exact syntax).  In
>>>>>>> fact,
>>>>>>>>>> you
>>>>>>>>>>> could even make "Geometric_Element" a keyword.  And if the user
>>>>>>> does
>>>>>>>>>> not
>>>>>>>>>>> define it, then you default to first order lagrange continuous
>>>>>>> element
>>>>>>>>>>> (which is what you do now).
>>>>>>>>>>>
>>>>>>>>>>> I think you have everything already to do this.  The
>>>>>>> Geometric_Element
>>>>>>>>>> is
>>>>>>>>>>> just a special case of a regular finite element, except the element
>>>>>>> is
>>>>>>>>>>> *always* the reference element.  You should be able to reuse all of
>>>>>>>>>> your
>>>>>>>>>>> FFC code generation.  Of course, computing higher order derivatives
>>>>>>> of
>>>>>>>>>>> basis functions will be more complicated (i.e. the chain rule will
>>>>>>>>>> give
>>>>>>>>>>> more terms b/c of the nonlinear map), but I assume you are doing
>>>>>>> most
>>>>>>>>>> of
>>>>>>>>>>> that symbolically within Python.
>>>>>>>>>>>
>>>>>>>>>>> 3. One last thing I learned in doing this was I also separated the
>>>>>>>>>>> transformations of the basis functions.  So the order of operations
>>>>>>> in
>>>>>>>>>> the
>>>>>>>>>>> assembly loop is:
>>>>>>>>>>>
>>>>>>>>>>> (a) get the index for the current element.
>>>>>>>>>>> (b) compute geometric transformations.
>>>>>>>>>>> (c) compute basis function transformations.
>>>>>>>>>>> (d) compute local element tensors.
>>>>>>>>>>> (e) insert into global matrix.
>>>>>>>>>>>
>>>>>>>>>>> At least for my setup, this gave a nice separation.  In fact, the
>>>>>>>>>> tabulate
>>>>>>>>>>> tensor routine is almost trivial since all the transformations are
>>>>>>>>>> done.
>>>>>>>>>>>
>>>>>>>>>>> I hope this info is helpful.  I am not trying to imply that my
>>>>>>> setup
>>>>>>>>>> is
>>>>>>>>>>> better; but maybe this gives a way for how to proceed.
>>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> Thanks. It's very helpful to hear your experience.
>>>>>>>>>>
>>>>>>>>>> I'm wondering if it's all as simple as adding the geometric info to
>>>>>>>>>> ufl::cell. It makes sense that a cell can fully describe its own
>>>>>>>>>> geometry.  We could have:
>>>>>>>>>>
>>>>>>>>>>  element = FiniteElement("triangle", 1, 1)
>>>>>>>>>>
>>>>>>>>>> which would be a general case in which a ufc::cell would be queried
>>>>>>> from
>>>>>>>>>> inside tabulate_tensor to provide the Jacobian, etc. For efficiency,
>>>>>>> we
>>>>>>>>>> could also have
>>>>>>>>>>
>>>>>>>>>>  element = FiniteElement("affine triangle", 1, 1)
>>>>>>>>>>
>>>>>>>>>> which would assume an affine map and therefore not query ufc::cell
>>>>>>> for
>>>>>>>>>> the transformation data.
>>>>>>>>>>
>>>>>>>>>> Garth
>>>>>>>>>>
>>>>>>>>>>> p.s. is it possible yet to have finite element spaces defined on a
>>>>>>>>>>> subdomain of codimension 1 simultaneously with other codimension 0
>>>>>>>>>> spaces
>>>>>>>>>>> and be able to compute bilinear forms involving mixed elements?  I
>>>>>>>>>> noticed
>>>>>>>>>>> you have some support for integrating on subdomains, but does that
>>>>>>>>>> include
>>>>>>>>>>> this???
>>>>>>>>>>>
>>>>>>>>>>> - Shawn
>>>>>>>>>>>
>>>>>>>>>>> On Wed, March 23, 2011 7:27 pm, Anders Logg wrote:
>>>>>>>>>>>> On Wed, Mar 23, 2011 at 11:24:37PM +0000, Garth N. Wells wrote:
>>>>>>>>>>>>> We've heard about this work, but never seen it ;).
>>>>>>>>>>>>> Can you post it somewhere? It could provide a discussion basis on
>>>>>>>>>> how
>>>>>>>>>>> to
>>>>>>>>>>>>> move this along.
>>>>>>>>>>>>> I think that we should be careful in asking FFC to do too much.
>>>>>>> It
>>>>>>>>>> may
>>>>>>>>>>> better to compute J, det(J), etc, outside tabulate_tensor(...)
>>>>>>>>>>>>> functions.
>>>>>>>>>>>>
>>>>>>>>>>>> Aren't those computed using codesnippets just pasted in by FFC
>>>>>>> into
>>>>>>>>>> the
>>>>>>>>>>> code anyway? Perhaps those codesnippets can be expanded with higher
>>>>>>>>>> order
>>>>>>>>>>> codesnippets?
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>>> Garth
>>>>>>>>>>>>> On 23/03/11 21:05, Peter Brune wrote:
>>>>>>>>>>>>>> I saw this thread and had already started running the tests I
>>>>>>> had
>>>>>>>>>> for
>>>>>>>>>>> this with the latest FFC to see if anything had changed recently.
>>>>>>> I
>>>>>>>>>> never
>>>>>>>>>>> made isoparametry through UFL public because I could never get
>>>>>>>>>>>>> the
>>>>>>>>>>>>>> efficiency to be reasonable.  The situation appears to have
>>>>>>>>>> improved
>>>>>>>>>>> a
>>>>>>>>>>>>>> little bit from a year ago, but not much.  With the latest FFC,
>>>>>>>>>> it's
>>>>>>>>>>> about 10x slower (in 2D) to assemble Poisson with P1 for the
>>>>>>>>>>>>> test/trial
>>>>>>>>>>>>>> space and P2 for the coordinates than the same problem with
>>>>>>> affine
>>>>>>>>>>> coordinates.  It's even worse if you include things like
>>>>>>>>>>>>>> Piola-transformed surface normals into the mix.  When I was
>>>>>>> working
>>>>>>>>>>> on
>>>>>>>>>>>>>> this I only assembled a very small fraction of the elements
>>>>>>> (around
>>>>>>>>>> a
>>>>>>>>>>> cylinder in a flow) with the parametric Jacobian, so it worked OK
>>>>>>> for
>>>>>>>>>>> small problems.
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> I believe that it would do much, much better if the optimizing
>>>>>>>>>>> quadrature compiler in FFC supported fractions.  This is necessary
>>>>>>>>>> because
>>>>>>>>>>> you need to apply the inverse isoparametric Jacobian, which
>>>>>>> includes 1
>>>>>>>>>> /
>>>>>>>>>>> |J|, to any basis function derivatives in the form.
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> For reference, here are the assembly times for the iso(super, I
>>>>>>>>>>>>> suppose
>>>>>>>>>>>>>> is more accurate)-parametric, optimized, and non-optimized
>>>>>>> Poisson
>>>>>>>>>>> problems in 2D on a 256x256 square (with the parametric coordinates
>>>>>>>>>>>>> just
>>>>>>>>>>>>>> being the regular coordinates passed through).
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> isoparametric assembly      |        3.2017      3.2017     1
>>>>>>>>>>> optimized assembly          |       0.34515     0.34515     1
>>>>>>> regular
>>>>>>>>>>> assembly            |       0.36524     0.36524     1
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> I got really deep in FFC trying to make this work with no
>>>>>>> success,
>>>>>>>>>>> but
>>>>>>>>>>>>>> this was before the rewrite.  I'd be willing to declare victory
>>>>>>> on
>>>>>>>>>>>>> this
>>>>>>>>>>>>>> one and submit my code if someone else were willing to make it
>>>>>>> fast
>>>>>>>>>>> enough to use.  There's also the issue of how exactly to extend the
>>>>>>>>>>> interface to support this in an elegant fashion.  Right now I just
>>>>>>>>>>>>> call
>>>>>>>>>>>>>> a function that takes a form and a parametric Jacobian, runs a
>>>>>>>>>>> Transformer on the form, and spits out a new form.
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> - Peter
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> On Wed, Mar 23, 2011 at 3:30 PM, Anders Logg <logg@xxxxxxxxx
>>>>>>>>>>>>>> <mailto:logg@xxxxxxxxx>> wrote:
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>     Peter Brune claims to have solved this by a small addition
>>>>>>> to
>>>>>>>>>> the
>>>>>>>>>>>>> form
>>>>>>>>>>>>>>     language that automatically expresses the curved elements as
>>>>>>> a
>>>>>>>>>>>>> mapping
>>>>>>>>>>>>>>     and expands appropriately (and invisible to the user) those
>>>>>>>>>>>>> mappings
>>>>>>>>>>>>>>     to yield a form that may then be assembled. The higher order
>>>>>>>>>>>>> geometry
>>>>>>>>>>>>>>     is then expressed as a vector-field on the mesh.
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>     Perhaps Peter can be pushed to polish up on his code and
>>>>>>> submit
>>>>>>>>>>>>> it.
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>     On Wed, Mar 23, 2011 at 07:46:57PM +0000, Garth N. Wells
>>>>>>> wrote:
>>>>>>>>>>>>>>     > We haven't really looked at this. It was discussed a while
>>>>>>>>>>> back,
>>>>>>>>>>>>>>     but no
>>>>>>>>>>>>>>     > one has committed much time to it. We struggled to settle
>>>>>>> on
>>>>>>>>>> an
>>>>>>>>>>> appropriate abstraction to push on with.
>>>>>>>>>>>>>>     >
>>>>>>>>>>>>>>     > Garth
>>>>>>>>>>>>>>     >
>>>>>>>>>>>>>>     > On 23/03/11 18:40, Douglas Arnold wrote:
>>>>>>>>>>>>>>     > > What is the status of curved (e.g., isoparametric)
>>>>>>> elements
>>>>>>>>>>> in
>>>>>>>>>>>>>>     dolfin?
>>>>>>>>>>>>>>     > > I gather they are not implemented in the main branch.
>>>>>>> Has
>>>>>>>>>>>>> anyone
>>>>>>>>>>>>>>     > > done anything with this can be used?  Is there any
>>>>>>> example
>>>>>>>>>>>>> code?
>>>>>>>>>>>>>>     > > (For example, if you want to
>>>>>>>>>>>>>>     > > solve the Poisson problem in a disc and get better than
>>>>>>> 2nd
>>>>>>>>>>>>> order
>>>>>>>>>>>>>>     > > convergence, you need to do better than polygonal
>>>>>>>>>>>>> approximation of
>>>>>>>>>>>>>>     > > the disc.)
>>>>>>>>>>>>>>     > >
>>>>>>>>>>>>>>     > >  -- Doug
>>>>>>>>>>>>>>     > >
>>>>>>>>>>>>>>     > > _______________________________________________
>>>>>>>>>>>>>>     > > Mailing list: https://launchpad.net/~dolfin
>>>>>>>>>>>>>>     <https://launchpad.net/%7Edolfin>
>>>>>>>>>>>>>>     > > Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
>>>>>>>>>>>>>>     <mailto:dolfin@xxxxxxxxxxxxxxxxxxx>
>>>>>>>>>>>>>>     > > Unsubscribe : https://launchpad.net/~dolfin
>>>>>>>>>>>>>>     <https://launchpad.net/%7Edolfin>
>>>>>>>>>>>>>>     > > More help   : https://help.launchpad.net/ListHelp
>>>>>>>>>>>>>>     >
>>>>>>>>>>>>>>     > _______________________________________________
>>>>>>>>>>>>>>     > Mailing list: https://launchpad.net/~dolfin
>>>>>>>>>>>>>>     <https://launchpad.net/%7Edolfin>
>>>>>>>>>>>>>>     > Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
>>>>>>>>>>>>>>     <mailto:dolfin@xxxxxxxxxxxxxxxxxxx>
>>>>>>>>>>>>>>     > Unsubscribe : https://launchpad.net/~dolfin
>>>>>>>>>>>>>>     <https://launchpad.net/%7Edolfin>
>>>>>>>>>>>>>>     > More help   : https://help.launchpad.net/ListHelp
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>> _______________________________________________
>>>>>>>>>>>> Mailing list: https://launchpad.net/~dolfin
>>>>>>>>>>>> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
>>>>>>>>>>>> Unsubscribe : https://launchpad.net/~dolfin
>>>>>>>>>>>> More help   : https://help.launchpad.net/ListHelp
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>
>>>>>>
>>>>>
>>>>>