Incomplete quadratures

[GoogleCodeExporter]

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

Original issue reported on code.google.com by freevryh...@gmail.com on 28 Apr 2010 at 3:02

[GoogleCodeExporter]

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.

Original comment by robert.c...@gmail.com on 28 Apr 2010 at 7:09

[GoogleCodeExporter]

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/

Original comment by robert.c...@gmail.com on 28 Apr 2010 at 3:07

[GoogleCodeExporter]

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.

Original comment by freevryh...@gmail.com on 29 Apr 2010 at 7:09

Attachments:

• 0002-Committer-freevryheid-freevryheid-gmail.com.patch

[GoogleCodeExporter]

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.

Original comment by robert.c...@gmail.com on 29 Apr 2010 at 8:06

[GoogleCodeExporter]

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?  

Original comment by freevryh...@gmail.com on 29 Apr 2010 at 12:28

[GoogleCodeExporter]

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?

Original comment by robert.c...@gmail.com on 29 Apr 2010 at 12:37

[GoogleCodeExporter]

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

Original comment by freevryh...@gmail.com on 29 Apr 2010 at 12:41

[GoogleCodeExporter]

OK, looking forward to it :)

Original comment by robert.c...@gmail.com on 29 Apr 2010 at 12:42

[GoogleCodeExporter]

Convert quadrature coordinates and weights to formulas when possible-

it's now a thing of beauty :)

Next stop is testing of the quadratures.

Original comment by freevryh...@gmail.com on 30 Apr 2010 at 1:06

Attachments:

• 0002-Convert-quadrature-coords-and-weights-to-formulas-wh.patch

[GoogleCodeExporter]

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

Original comment by logan.so...@gmail.com on 30 Apr 2010 at 1:17

[GoogleCodeExporter]

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!

Original comment by robert.c...@gmail.com on 30 Apr 2010 at 12:31

[GoogleCodeExporter]

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

Original comment by robert.c...@gmail.com on 3 May 2010 at 1:38

[GoogleCodeExporter]

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.

Original comment by robert.c...@gmail.com on 4 May 2010 at 10:24

[GoogleCodeExporter]

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).

Original comment by robert.c...@gmail.com on 6 May 2010 at 11:39

[GoogleCodeExporter]

This can be reopened if we need higher orders later.

Original comment by freevryh...@gmail.com on 18 May 2010 at 3:29

• Changed state: Fixed

[GoogleCodeExporter]

Migrated to http://github.com/sfepy/sfepy/issues/127

Original comment by robert.c...@gmail.com on 30 Jan 2012 at 10:26

• Changed state: Migrated