sfepy/sfepy — Issue #127: Incomplete quadratures
Author: rc
Created: 2012-01-30 10:26:52+00:00
State: open

The quadratures are incomplete and some may not have not the order they
claim. These need to be reviewed.

Migrated from http://code.google.com/p/sfepy/issues/detail?id=125

earlier comments

robert.cimrman said, at 2010-04-28T07:09:50.000Z:

Oh yes, thanks for creating this issue.

IMHO we should refactor sfepy/fem/quadratures.py - I would store the quadrature
weights, points, etc. in "tables" using dictionaries (just like we store element
geometries in sfepy/fem/geometry_element.py) and put all the functionality in a
single class Quadrature.

Not sure yet what kind of keys/hierarchy of dicts would work well.

robert.cimrman said, at 2010-04-28T15:07:09.000Z:

I have refactored integrals and quadratures, see [1]. Now it remains to write some
tests into tests/test_quadratures.py, and add more quadratures.

[1] http://github.com/rc/sfepy/

freevryheid said, at 2010-04-29T07:09:29.000Z:

Wow, much simpler now!!

Validated existing, fixed and added to orders of quad coords + weights. See attached
patch.

Followed the docs but not sure I did this correctly. Need sleep :)

Referenced "The Finite Element Method Displayed" by Dhatt and Touzat - knew it would
come in handy one day :)

These still need to be tested.

robert.cimrman said, at 2010-04-29T08:06:26.000Z:

Nice! I just fixed the commit message and author name - see the github repo - is it ok?

BTW. I suggest you put a .gitconfig into your home with

[user]
email = freevryheid...
name = Andre Smit

BTW2. The references could be mentioned in the module docstring.

freevryheid said, at 2010-04-29T12:28:02.000Z:

Thinking about it I realized that it wouldn't make a significant difference if the
coords and weights are input as formulae e.g. 1.0/np.sqrt(6.0) rather than the result
i.e. 0.40824829046386307. I've used the latter but prefer the former. I'm thinking
about updating. Any thoughts?

Also, is there any benefit to including higher orders?

robert.cimrman said, at 2010-04-29T12:37:04.000Z:

Higher orders might be good to have in future. It would be desirable to be able to
set the 'automatic' order, that would be determined by each term, given polynomial
orders of its arguments.

As for using formulas, they might provide more insight, yes, and incur no runtime
penalty, so why not. On the other hand, there are coordinates where no such nice
formulas exist. So do as you wish.

Will you create a new patch w.r.t. the current HEAD at github, or you prefer amending
your patch?

freevryheid said, at 2010-04-29T12:41:18.000Z:

Let's use formulae if we can - I'll create a new patch this evening. Thanks for git
help.

robert.cimrman said, at 2010-04-29T12:42:40.000Z:

OK, looking forward to it :)

freevryheid said, at 2010-04-30T01:06:31.000Z:

Convert quadrature coordinates and weights to formulas when possible-

it's now a thing of beauty :)

Next stop is testing of the quadratures.

logan.sorenson said, at 2010-04-30T01:17:31.000Z:

+1 I like the formula approach, if just for improving the appearance of the code! :)

robert.cimrman said, at 2010-04-30T12:31:36.000Z:

It's at [1] now. I have just removed the numpy import, as that is achieved by the
'from sfepy.base.base import *' statement

Thanks!

robert.cimrman said, at 2010-05-03T13:38:05.000Z:

So a basic test of quadratures is now at [1]. Next steps are:

- write QuadraturePoints.integrate() or Integral.integrate(), that would compute the
integral for a given function of space, passed as argument, over a domain given by
the geometry key. Probably the latter (Integral.integrate()) would be better, as
Integral knows about the geometry, but QuadraturePoints could be enhanced to do so too.

- make up some functions of given polynomial order, and add test_quadratures() to
test the quadratures using the above function

robert.cimrman said, at 2010-05-04T10:24:30.000Z:

So it's at [1] - all quadratures pass the test except '2_4', order 5. Could you
verify/correct that? Then we can close this issue.

robert.cimrman said, at 2010-05-06T11:39:19.000Z:

The following composite Gauss quadrature works:

5 : _QP([[-nm.sqrt(3.0/5.0), -nm.sqrt(3.0/5.0), (5.0/9.0) * (5.0/9.0)],
[ 0.0 , -nm.sqrt(3.0/5.0), (8.0/9.0) * (5.0/9.0)],
[ nm.sqrt(3.0/5.0), -nm.sqrt(3.0/5.0), (5.0/9.0) * (5.0/9.0)],
[-nm.sqrt(3.0/5.0), 0.0, (5.0/9.0) * (8.0/9.0)],
[ 0.0 , 0.0, (8.0/9.0) * (8.0/9.0)],
[ nm.sqrt(3.0/5.0), 0.0, (5.0/9.0) * (8.0/9.0)],
[-nm.sqrt(3.0/5.0), nm.sqrt(3.0/5.0), (5.0/9.0) * (5.0/9.0)],
[ 0.0 , nm.sqrt(3.0/5.0), (8.0/9.0) * (5.0/9.0)],
[ nm.sqrt(3.0/5.0), nm.sqrt(3.0/5.0), (5.0/9.0) * (5.0/9.0)]],
bounds=(-1.0, 1.0)),

but it needs 9 points. Do you think you 8-point scheme could be fixed? (I googled,
but found nothing).

freevryheid said, at 2010-05-18T03:29:00.000Z:

This can be reopened if we need higher orders later.