pombreda / sfepy

Automatically exported from code.google.com/p/sfepy
0 stars 0 forks source link

Incomplete quadratures #125

Open GoogleCodeExporter opened 9 years ago

GoogleCodeExporter commented 9 years ago
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 commented 9 years ago
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 commented 9 years ago
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 commented 9 years ago
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:

GoogleCodeExporter commented 9 years ago
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 commented 9 years ago
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 commented 9 years ago
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 commented 9 years ago
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 commented 9 years ago
OK, looking forward to it :)

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

GoogleCodeExporter commented 9 years ago
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:

GoogleCodeExporter commented 9 years ago
+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 commented 9 years ago
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 commented 9 years ago
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 commented 9 years ago
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 commented 9 years ago
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 commented 9 years ago
This can be reopened if we need higher orders later.

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

GoogleCodeExporter commented 9 years ago
Migrated to http://github.com/sfepy/sfepy/issues/127

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